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Abstract — A  promising  approach  to  detecting  roadside  bombs  attached  to  com¬ 
mand  wires  is  the  electromagnetic  sensing  and  identification  of  the  wires.  The  lowest 
five  resonant  frequencies  of  the  wires,  along  with  the  widths  of  the  resonances,  can 
serve  as  a  “fingerprint”  for  finding  the  wires.  A  first  large  step  toward  exploiting  this 
fingerprint  is  to  calculate  the  resonances  and  their  widths  for  a  straight  wire  on  a  flat 
interface  between  a  homogeneous  earth  and  air.  The  calculation  of  resonances  re¬ 
quires  extending  the  theory  of  the  linear  antenna  to  deal  with  a  wire  on  the  interface 
between  two  dielectric  media,  which  we  accomplish  here.  Complex-valued  resonant 
frequencies  are  defined  as  those  for  which  a  certain  homogeneous  integral  equation 
for  the  current  in  the  strip  on  the  interface  has  non-trivial  solutions.  By  applying 
a  Galerkin  procedure  we  obtain  approximate  numerical  solutions  for  the  resonant 
frequencies  and  their  widths.  A  table  of  resonances  is  given  for  the  case  of  a  relative 
dielectric  constant  of  the  earth  equal  to  4  and  for  three  values  of  the  ratio  of  wire 
radius  a  to  wire  length  h.  MATLAB  computer  programs  for  determining  resonant 
frequencies  and  widths  for  other  parameter  values  are  included. 
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I.  INTRODUCTION  AND  SUMMARY  OF  RESULTS 


The  need  to  detect  and  identify  command  wires  that  control  roadside  bombs  raises  the 
problem  of  determining  an  electromagnetic  signature  of  command  wires  that  run  along 
the  surface  of  the  earth.  We  study  such  a  wire  (or  wire  pair)  in  terms  of  the  theory  of 
linear  antennas.  The  theory  of  linear  antennas  in  a  single  homogeneous  medium  was  well 
developed  over  hfty  years  ago  in  terms  of  integral  equations  [1];  however,  the  problem  of 
characterizing  command  wires  calls  for  extending  this  theory  to  deal  with  two  media  at 
a  planar  interface,  and  the  extra  medium  raises  a  whole  new  set  of  problems,  requiring  a 
substantial  development  of  technique. 

In  terms  of  the  Maxwell  equations  that  dehne  electromagnetism,  here  we  report  on 
the  needed  technique  and  its  application  to  the  determination  of  resonant  frequencies  of 
a  straight  thin  wire  (or  wire  pair)  on  a  flat  interface  between  earth  and  air  or,  more  gen¬ 
erally,  between  two  homogeneous  lossless  media,  each  characterized  by  a  given  dielectric 
constant.  In  case  of  a  thin  wire  in  a  homogeneous  medium,  the  shape  of  the  cross  section 
is  unimportant.  For  example,  with  perfect  conductors,  a  wire  consisting  of  a  flat  strip  of 
width  4a  and  negligible  thickness  is  electromagnetically  essentially  the  same  as  a  circular 
wire  of  radius  a.  This  insensitivity  to  cross-sectional  shape  holds  also  for  a  thin  wire  or 
wire  pair  on  the  boundary  between  two  different  dielectric  media.  Therefore  one  is  free  to 
pick  a  shape  of  the  cross  section,  and  the  analytically  convenient  shape  is  a  flat,  perfectly 
conducting  strip  of  length  2h,  width  4a,  and  zero  thickness.  We  assume  the  strip  is  located 
in  a  plane  interface,  thought  of  as  horizontal,  between  a  medium  below  (earth)  with  a  di¬ 
electric  constant  ki  and  a  medium  above  (air)  with  a  dielectric  constant  k2-  The  problem 
is  to  hnd  the  hrst  hve  resonant  frequencies  for  electromagnetic  radiation  scattered  by  this 
thin  strip  as  functions  of  the  parameters  a,  h,  /ci,  and  /c2;  the  width  of  each  resonance  is 
also  to  be  determined.  For  a  single  homogeneous  medium,  our  understanding  both  of  the 
linear  antenna  and  the  related  scattering  problem  comes  mostly  from  the  Hallen  integral 
equation.  From  this  integral  equation,  in  which  the  vector  potential  plays  a  central  role,  it 
follows  that  the  current  on  such  an  antenna  is  approximately  sinusoidal,  especially  when  the 
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antenna  is  thin.  Fnrthermore,  the  wave  nnmber  for  this  approximately  sinnsoidal  cnrrent  is 
that  of  the  snrronnding  medinm.  In  contrast,  for  the  present  case  of  a  linear  antenna  on  an 
interface,  there  is  no  nsefnl  definition  of  a  vector  potential,  which  complicates  formnlating 
the  problem.  When  the  wire  involves  two  media  with  distinct  propagation  constants,  the 
qnestions  arise: 

1.  Is  the  cnrrent  along  the  wire  approximately  sinnsoidal;  and 

2.  if  so,  what  is  its  wave  nnmber? 

Based  on  the  known  expression  for  the  electric  field  emanating  from  a  point  cnrrent  element 
located  on  and  parallel  with  the  interface  [2] ,  integral  eqnations  of  both  the  Pocklington  type 
and  the  Hallen  type  can  be  formnlated,  bnt  with  a  more  complex  kernel  that  introdnces 
qnalitatively  new  featnres,  associated  physically  with  the  complexity  of  paths  by  which 
energy  can  propagate  near  an  interface.  As  described  below,  it  tnrns  ont  that  the  cnrrent 
near  resonance  is  approximately  sinnsoidal,  bnt,  even  when  both  media  are  lossless,  with  a 
comp  lex- valned  propagation  constant. 

Even  for  a  single  medinm  the  Hallen  and  the  Pocklington  eqnations  for  antenna  problems 
are  notorionsly  difficnlt  to  solve  either  analytically  or  nnmerically,  and  indeed  have  been 
solved  accnrately  only  for  a  few  special  cases,  snch  as  the  center-driven  antenna  [1] .  Recently, 
however,  we  defined  resonances  in  terms  of  special  freqnencies  in  the  complex  plane  at  which 
the  homogeneons  Pocklington  eqnation  has  non-zero  solntions.  We  then  determine  the  first 
few  resonances  approximately  by  nse  of  a  Galerkin  method,  withont  the  necessity  of  solving 
any  integral  eqnation.  We  checked  that  the  application  of  this  method  to  the  case  of  a  wire 
in  a  single  homogeneons  medinm  showed  agreement  with  known  resonances  given  in  Ref.  [1] 
(in  which  the  cnrrent  distribntion  along  the  wire  is  symmetric). 

We  proceed  as  follows: 

1.  In  Sec.  II,  taking  the  strip  to  rnn  along  the  x-axis,  we  start  with  the  expression  of 
the  electric  field  Ex{x,y)  along  the  strip  generated  by  a  cnrrent  element  on  the  strip, 
obtained  as  the  inversion  of  a  Fonrier  transform.  This  expression  lays  the  gronndwork 
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for  a  kernel  in  an  integral  equation,  once  one  finds  a  way  to  avoid  a  non-integrable 
singularity. 

2.  In  Sec.  Ill  we  find  the  complex-valued  zeros  in  the  Fourier  transform  of  Ex{x,y)  with 
respect  to  x,  denoted  Ex{(,y).  These  zeros  play  a  role  similar  to  the  zeros  in  the 
Fourier  transform  for  the  homogeneous  case,  and  in  the  same  way  they  allow  one 
to  obtain  an  integral  equation  of  the  Pocklington  type  for  the  current  in  the  strip, 
arranged  to  evade  non-integrable  singularities. 

3.  In  Sec.  IV  we  derive  a  suitably  symmetrized  integral  equation  of  the  Pocklington  type 
for  the  current  in  a  thin  strip  on  the  interface.  Important  implications  are:  (1)  that 
the  current  along  the  wire  near  resonance  is  indeed  sinusoidal,  and  (2)  we  have  the 
propagation  constant  of  this  approximately  sinusoidal  current. 

4.  In  Sec.  V  we  dehne  the  complex-valued  resonant  frequencies  as  the  frequencies  at 
which  the  homogeneous  integral  equation  for  the  current  in  the  strip  has  non-zero 
solutions.  Approximate  equations  for  the  resonant  frequencies,  suitable  for  numerical 
calculation,  are  obtained  by  use  of  a  Galerkin  procedure. 

5.  In  Sec.  VI,  the  equations  are  rearranged  to  facilitate  numerical  computation,  largely 
by  changing  integration  contours  to  avoid  undue  oscillation  in  the  integrands  of  various 
integrals. 

6.  Appendix  F  gives  the  MATLAB  programs  that  generate  numbers  from  the  rearranged 
equations  to  provide  numerical  values  for  the  hrst  hve  resonant  frequencies,  along 
with  the  width  of  each  resonance  (which  is  simply  related  to  the  imaginary  part  of  the 
complex- valued  resonant  frequency). 

As  an  example  of  resonant  frequencies,  here  are  the  results  from  the  analysis  given  below 
(computed  with  the  MATLAB  programs  in  Appendix  F)  for  the  case  of  earth  having  a 
dielectric  constant  4  times  that  of  free  space,  leading  to  /ci//c2  =  2.  The  hrst  hve  complex¬ 
valued  resonant  frequencies  for  three  cases  of  a/h,  the  ratio  of  wire  thickness  to  wire  length. 
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are  given  in  Table  1.  The  propagation  constant  in  air  is  k2  =  271  f  /c  where  /  is  the  frequency 
in  Hertz  and  c  ~  3  x  10®  m/s  is  the  speed  of  light.  (Note  that  an  even  resonance  number  n 
corresponds  to  anti- resonance  for  the  impedance  of  a  center-driven  antenna.) 

Figures  1,  2,  and  3  show  the  resonances  implied  by  Table  1  as  a  function  of  frequency. 


TABLE  1:  Complex  values  of  k2h  at  resonance  n 


n 

a//i  =  10  ^ 

a//i  =  10-5 

a/h  =  10-6 

1 

0.956824-i0.047401 

0.965436-i0.038153 

0.970820-i0.031885 

2 

1.934017-i0.079959 

1.947398-i0.063718 

1.955503-^0.052866 

3 

2.905922-i0. 115751 

2.925628-i0.091953 

2.937355-i0.076062 

4 

3.872673-Z0.153152 

3.900158-i0.121734 

3.916394-^0.100618 

5 

4.837435-iO.  190228 

4.873157-i0.151508 

4.894259-iO.  125246 

FIG.  1:  Resonant  scattering  from  wire;  6^=4,  a//i  =  10 
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II.  FORMULATION 

The  behavior  of  electromagnetic  waves  along  and  across  the  boundary  between  two  me¬ 
dia  with  different  properties  is  much  more  complicated  than  the  case  of  one  homogeneous 
medium.  Fortunately,  this  case  of  two  media  has  been  systematically  studied  [2],  The 
notations,  approach,  and  results  of  this  reference  are  to  be  used  here. 

For  dehniteness,  let  the  two  media  be  air  and  earth.  Following  the  notation  of  Ref.  [2] , 
let  the  planar  boundary  between  air  and  earth  be  in  the  a;|/-plane  (i.e.,  =  0)  and  the  z-axis 

points  in  the  direction  of  the  earth.  The  problem  of  the  strip  antenna  to  be  studied  is  shown 
schematically  in  Fig.  4.  The  length  of  the  antenna  is  2h.  Still  following  the  notation  of  Ref. 
[2],  let  ki  be  the  wave  number  of  region  1  (earth)  and  /c2  that  of  region  2  (air).  As  holds  for 
earth  and  air,  it  will  be  assumed  throughout  this  report  that 

h  >  k2.  (2.1) 
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Since  the  antenna  is  assumed  to  be  thin,  a  is  taken  to  be  small,  so  that 

a  h  and  kia  1,  (2.2) 

which,  of  course,  implies  k2a  -C  1.  For  present  purposes  both  earth  and  air  are  assumed 
lossless,  so  that  both  ki  and  k2  are  taken  to  be  real-valued  and  positive.  As  in  the  usual 
case  of  the  linear  antenna  in  a  homogeneous  medium,  once  the  lossless  case  is  understood, 
the  introduction  of  loss  in  the  medium  is  fairly  straightforward. 

The  formulation  of  the  problem  hinges  on  the  electric  held  generated  by  a  “horizontal” 
point  dipole  in  the  interface.  Consider  a  delta-function  current  at  the  origin. 


J  =  S(x)6(y)S{z)x-, 


(2,3) 


then,  the  x-component  of  the  electric  held  at  a  point  {x,  y,  0)  on  the  interface  is  given  by 
(see  Eq.  (5.4.13)  of  Ref.  [2]) 


Ex{x,y)  = 


cupo 


/  dye 


X^x+r^y)  7l(fei  +72(^1  -^^) 


'  —00  J  —CO 


MN 


(2.4) 
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where 


7i  = 


72  =  {ki  - 


(2.5) 


M  =  7i  +  72,  N  =  kl'j2  + 

In  deriving  Eq.  (2.4),  it  has  been  assumed  that  the  /i  for  both  region  1  and  region  2  is  given 
by  /io;  that  is,  we  do  not  deal  here  with  magnetizable  earth.  Note  that  Eq.  (2.4)  does  not 
change  when  the  regions  1  and  2  are  exchanged. 

The  advantage  of  studying  a  strip  antenna  of  negligible  thickness  can  be  seen  from 
Eq.  (2.4).  If  the  antenna  is  circular  for  example,  then  both  the  current  source  (2.3)  and 
the  observation  point  would  have  non-negligible  values  for  the  ^-coordinate.  In  this  case, 
the  corresponding  formula  depends  on  the  signs  of  these  2;  coordinates  and  is  much  more 
complicated  than  the  present  Eq.  (2.4). 


x  =  -h 

1 

1 

1 

1 

1 

1 

1 

Region  2:  air 

x=  h 

X 

1 

strip  antenna  of 

1 

negiigibie  thickness 

1 

1 

1 

1 

Region  1 :  earth 

I  z 


(a)  side  view 


x  =  -h 


x  = 


V 


y=  -2a 


X 


y=  2a 


I  y 

(b)  top  view 


FIG.  4:  Schematic  diagram  of  the  antenna  geometry. 
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III.  ZEROS  IN  THE  FOURIER  TRANSFORM  OF 


The  current  induced  in  the  strip  will  be  shown  to  satisfy  an  integral  equation  of  the 
Pocklington  type.  In  the  case  of  a  single  medium,  as  outlined  in  Appendix  A,  the  kernel 
for  the  integral  equation  is  derived  easily  once  one  notices  a  pair  of  zeros  in  the  Fourier 
transform  of  E^.  That  method  works  also  for  for  the  case  of  a  thin  strip  at  an  interface, 
except  that  one  has  to  do  some  work  to  find  the  zeros.  It  is  to  this  task  that  we  now  turn. 

From  Eq.  (2.4)  it  follows  that  the  Fourier  transform  from  a;  to  of  Ex{x,  y)  is  proportional 


to 


my)  -  r 

J  —oo 


e)+i2{ki 

MN 


This  E{^,y)  depends  also  on  ki  and  /c2,  which  in  turn  depend  on  frequency  through 


(3.1) 


ki  =  ^/Eiuj/c  and  /c2  =  ^/oiui/c, 


(3.2) 


where  c  is  the  speed  of  light  and  ei  and  62  are  the  relative  dielectric  constants  of  the  two 
media  (taking  region  2  to  be  air,  /c2  ~  1).  We  assume  lossless  media,  so  that  both  ei  and  62 
are  real-valued  and  positive.  In  marked  contrast  to  the  case  of  a  single  medium,  however, 
it  will  turn  out  that  the  zeros  of  E{^,  y)  will  be  complex- valued.  To  proceed,  we  re-write 
Eq.  (3.1)  as 

my)--eFm)+m.y).  (3.3) 


where 


Fi{ty)=  /  dye 


irjy 


and 


F2{i,y)=  /  dye 


kly2  +  kjyi 

1 


irjy 


(3.4) 


(3.5) 


7i  +  72 

Of  these  two  functions,  E2  is  the  simpler  one;  in  fact,  it  can  be  expressed  exactly  in  terms 
of  Bessel  functions: 

TT  1 


F2iC,y)  = 


kl-klyl 


,/]meH\^\y^/Wm)-^/WmHY\y^/Wm)  ■  (3.6) 


Unfortunately,  there  does  not  seem  to  be  a  corresponding  formula  for  Ei.  However,  such 
an  exact  evaluation  for  Ei  is  not  needed.  Since  the  thin  antenna  is  by  far  the  most  important 


case,  as  expressed  by  the  condition  /ca  1  of  (2.2),  there  is  a  corresponding  condition  for 
this  Fi.  Here  y  plays  the  role  of  a,  because  ||/|  <  4a.  Taking  y  to  be  positive  without  loss  of 
generality,  the  corresponding  condition  is 


kiy  <  1  and  |,^||/  -C  1. 


(3.7) 


An  absolute  value  sign  has  been  kept  on  ^  because,  as  to  be  seen  later,  ^  is  allowed  to  be 
complex. 

Before  evaluating  Fi  under  condition  (3.7),  it  is  instructive  to  write  down  the  correspond¬ 
ing  approximate  formula  for  the  F2  of  Eq.  (3.6).  This  is  immediately 


TT 


F2iC,y)  ~  2 


{kj  -  ^2)  \n{y^kf  -  )  -  (/ci  -  C^)  )  ,  1  o  ^ 

- - +7-ln2-  - 


(3.8) 

Furthermore,  from  the  series  expansion  for  the  Bessel  function,  the  first  neglected  term  is 
of  the  order  of  {CuYl  where  ^  is  the  larger  one  of  ki  and  |.^|.  This  point  turns  out 

to  be  of  importance:  what  it  means  is  that  approximation  (3.8)  is  highly  accurate  for  small 
values  of  y. 

The  series  expansion  for  this  F2{^,y)  is  of  the  form 


ln(Cl/),  1,  iCyfHCy),  (C2/)^  (C2/)^ln(C2/),  (Cl/)^  •••  (3-9) 


for  small  values  of  y.  As  proved  in  Appendix  B  the  corresponding  series  expansion  for  the 
Fi{^,y)  is  also  of  the  form  (3.9).  The  task  is  therefore  to  determine  the  first  two  terms, 
namely  those  of  orders  ln(,^|/)  and  1,  for  Fi{^,y)  from  Eq.  (3.4). 

As  a  first  step  in  this  direction,  let  Fi  y)  be  written  in  the  following  form: 


Fi{U)  = 


-1 


kl  -  kl 

The  denominator  here  has  a  zero  at 


-  (if  +  i|)(e  +  if 


(3.10) 


y +1) 


2 


^1^2 

kl  +  kl' 


(3.11) 


This  is  the  famous  Sommerfeld  pole  [3].  At  this  point,  the  numerator  on  the  right-hand  side 
of  Eq.  (3.10)  also  has  a  zero,  implying  the  well-known  fact  that  the  Sommerfeld  pole  is  on 
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the  second  sheet.  From  either  Eq.  (3.4)  or  Eq.  (3.10),  it  is  readily  seen  that  the  last  factor  in 
the  integrand  is  of  order  77“^  for  large  values  of  rj.  Therefore,  the  leading  term  of  order  ln{^y) 
for  small  values  of  y,  as  discussed  in  the  preceding  paragraph,  comes  from  the  integration 
over  large  values  of  rj.  It  is  thus  desirable  to  separate  out  the  part  of  the  integrand  that 
behaves  as  for  large  r];  in  the  rest  of  the  terms,  y  can  be  simply  put  to  zero.  One  way  to 
accomplish  this  is  to  rewrite  the  right-hand  side  of  Eq.  (3.10)  in  the  somewhat  complicated 
form  as  follows: 


Fi{U)  = 


^2  _  U2  U2  \  U2  / 

klk 


drj  e 


ir]y 


kl 


kl 


7]^ 


+ 


27.2 
2 


kl 


kl 


(3.12) 


kiki  -  [kl  +  kiw + \ ^ki  -e-v^  Vki  -e-v^ 

This  way  of  rewriting  Eq.  (3.10)  has  several  advantages.  First,  each  of  the  four  terms  on  the 
right-hand  side  contains  either  \/kl  —  —  7f  or  \/kl  —  —  rf,  but  not  both.  Secondly,  the 

first  two  integrations  lead  to  Bessel  functions.  Thirdly,  the  quantity  within  the  parentheses 
is  zero  at  the  Sommerfeld  pole  given  by  (3.11).  Moreover,  the  integrand  with  this  quantity 
as  a  factor  decreases  as  ri~^  for  large  rj,  and  hence  y  can  be  put  to  zero  in  this  integral, 
simplifying  its  evaluation  greatly.  These  considerations  lead  to 
1 


Fiity) 


TT 


1-2  I  l2  )  U2  _  U2 


+  F, 


30 


(3.13) 


where 

-1  r,  m  (  kl  kl 

kl  —  kl  klkl  —  {kl  +  kl){^‘^  -f-  jf)  ^ —  ^2  _  ^2 

(3.14) 

is  independent  of  y.  Approximation  (3.13)  is  accurate  enough  to  give  for  Fi  the  two  leading 
terms  of  orders  ln(^|/)  and  1. 

Since  only  these  two  leading  terms  are  to  be  calculated,  the  Bessel  functions  on  the 
right-hand  side  of  Eq.  (3.13)  can  be  expanded  to  give 

+  -fso  , 

(3.15) 


Fi{U) 


kl  kl 


Ti  +  2i 


kl  ln{y^/^FF^)  -  kl  \n{y^/^Fr^ 

7^1  7^2 


-|-  7  —  In  2 
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while  a  direct  evaluation  of  Eq.  (3.14)  gives 


E?n  —  ■ 


k'i-kj  ^(k:f  +  kl)e  -  k\kl 


kf + + k'i)e  -  njki  _  fc|  +  ^(kkj  +  ki)e  -  kinj 
ki-^(ki  +  i4)e-kiki  "  ki  -  ^ik\  +  q)e  -  kie. 

Equation  (3.16)  has  been  written  down  for 

^  f-2  ^  >2 


(3.16) 


/cf  +  /c| 


<  r  <  K. 


(3.17) 


and  analytic  continuation  can  be  used  to  ascertain  the  values  of  E30  outside  this  range  of 
In  summary,  the  substitution  of  (3.15),  (3.16)  and  (3.8)  into  (3.3)  gives  explicitly 


,  .  (/c?  -  C^)  \n{y^kl  -  )  -  (/ci  -  C^)  \ri{y^kl  -  ) 


E  ~ - \r  i 

2 


[  kl-kl 

I  •  -^2)  _  /g2  ln(|/y^/c2  _^2) 

l2  I  l2  ■)  o  ^  U2  _  U2 


+  7  —  In  2 - 

^  2 


+  7  —  In  2 


ikfkl  1 

~W^27WT^W^WM 

^  k?  +  +  ki)e  -  kiki  _  +  ^(ki  +  ki)e  -  kiki 


(3.18) 


This  is  the  desired  approximate  expression  from  which  to  determine  the  zero  of  this  E  as  a 
function  of  ^  when  ki,  /c2  and  y  are  given  subject  to  the  conditions  (3.7). 

As  a  zeroth  approximation,  only  the  leading  term  of  order  ln(,^|/)  is  kept  while  all  terms  of 
order  1  are  neglected.  In  this  approximation,  the  right-hand  side  of  (3.18)  is  much  simplified: 


E  ^  i  [  1  — 


kl  +  k^ 


lii(?!/). 


(3.19) 


Thus,  in  this  rough  approximation,  the  zero  of  E,  called  ko,  is  given  simply  by 

^2  I  u2 
u2  W  ^2 


(3.20) 


What  Eq.  (3.20)  means  physically  is  as  follows.  Since  this  location  ko  of  the  zero  of  E 


is  real,  the  current  on  a  linear  antenna  at  the  planar  boundary  of  two  lossless  dielectric 
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materials  is  approximately  sinusoidal  with  a  wave  number  given  by  a/ {kf  +  /c|)/2.  This 
answers  the  two  questions  posed  in  the  Introduction. 

It  must  be  remembered  that  (3.20)  is  only  a  rough  approximation.  It  is  the  purpose  of 
the  remainder  of  this  section  to  obtain  a  better  approximation  by  taking  into  account  not 
only  the  leading  term  of  order  ln{^y)  but  also  the  next-to-leading  term  of  order  1.  Since  the 
right-hand  side  of  (3.18)  is  complex  in  the  sense  of  having  both  a  real  part  and  an  imaginary 
part,  setting  E  to  zero  can  be  expected  to  lead  to  a  solution  ^  that  is  complex.  Such  a  zero 
off  the  real  axis  leads  to  many  new  phenomena. 

The  zero-th  order  approximation  (3.20)  implies  that,  to  this  approximation. 


On  the  basis  of  (3.21),  define  the  dimensionless  quantity 


ki  -  k. 


2  g*7r/2 


y  =  y 


_  U2 


and  (3.18)  takes  the  form 

,ki+ki-2e 


E 


\n.y  +  ^ 


ln2-  — 
2 


+  i 


kj  +  kl 

(k't-e)Ha  ln(a  ) 


2i(‘‘ 


k(  +  /C2 

^2l2 


1.2  _  l2 

kl  \n{a^/kl  -^2)  _  j.2  \^(^oikJkl  - 

/il  /i2 


kl  -  kl  ^[kl  -h  kl)e  -  kjkl 


X 


+  ywTfcpe  -  +  q)e  -  kjki 


\/*¥^ 


where 


a  = 


(3.21) 


(3.22) 


(3.23) 


(3.24) 


U2  _  k2  ' 

In  (3.23),  the  small  parameter  y  appears  only  in  the  very  first  term.  In  the  notation  of  (3.9), 
this  first  term  is  the  only  one  of  order  ln{^y),  while  everything  else  is  of  order  1. 
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It  is  curious  to  note  that  the  a  of  Eq.  (3.24)  is  infinity  when  ki  =  /c2,  i.e.,  when  the  two 
media  have  the  same  electromagnetic  properties.  In  this  case,  the  zero  of  E  is  given  exactly 
by  C  =  (=  ^2),  as  already  mentioned  in  Sec.  III.  Nevertheless,  the  case  where  ki  and  /c2 

are  close  to  each  other  requires  a  careful  treatment.  Since  this  case  is  not  of  great  practical 
importance,  it  will  be  dealt  with  in  Appendix  C 

The  approximation  (3.23)  can  be  used  to  get  a  better  approximation  to  the  zero  of  E 
than  that  given  by  (3.20).  A  major  reason  to  seek  a  better  approximation  is  to  find  the 
imaginary  part  of  this  zero,  including  whether  it  remains  zero.  For  this  purpose,  define 


^(0)  =  Eiy) 

'\e={kl+kl)/2 

When  (3.24)  is  used  for  E{^,y),  then  this  is  independent  of  y. 
It  follows  from  (3.21)  that 


(3.25) 


\n(a\/kl  — 

Accordingly,  this  E^^^  is  given  by 


=  0  and  \n{a\/k2  —  ) 


m 

e=(kl+kl)/2  ~  2  ■ 


(3.26) 


^(0) 


m 

~4 


T 

2 

2l2 


kl 


kik 


kl 

1 


m 


kl  2 


j  /c^a/2  +  ^Jk\  +  kl  ^  in 


kl  -  kl  ^(kf  +  kl)/2 

Note  that  this  E^^'>  has  both  real  and  imaginary  parts. 

In  terms  of  this  E^^\  the  zero  of  E{^,y)  is  given  approximately  by 


(3.27) 


kl- 


kl  +  kl 


(3.28) 


\n{yy^{kl  -  kl)/2  )+  7  -  In  2  -  in {2 _ 

In  the  denominator,  the  terms  7— In  2  and  — i7r/2  may  or  may  not  be  kept.  Roughly  speaking, 
it  is  the  real  part  of  E^^'^  that  determines  the  imaginary  part  of  the  zero  of  E{y). 

Combining  the  real  and  imaginary  parts  separately,  the  E^^'>  as  given  in  (3.9)  can  be 
rewritten  as 


9  h.2  _  h2 


kl  +  k 


klk 


2u2 
2 


\/  {k\  +  kl)/2 


—  i 


klk 


2u2 
2 


kl\/2  +  \/kf  kl 
kl  -  kl  ^{kf  +  kl)/2  klV2  + 


In 


(3.29) 
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TABLE  2:  Sample  values  of  the  as  a  function  of  k2a  and  k\/k2- 


k2a 

ki/k2  =  1.8 

ki/k2  =  2.0 

ki/k2  =  2.2 

10-^ 

2.18318 -hiO.12600 

2.59013 -h  i0.16384 

3.04044  -h  i0.20489 

10-5 

2.16737 -hiO.10080 

2.56777 -h  i0.13092 

3.01066 -h  i0.16348 

10-6 

2.15777 -hiO.08392 

2.55413 -h  i0.10890 

2.99250  -h  i0.13586 

For  the  real  part,  since 


1 


it  follows  that 


< 


2/c^/c|  kl  +  k. 


v'(*l  +  ^  *^1  +  *=2 


< 


and,  hence,  on  the  basis  of  (3.29), 


Re  >  0 


(3.30) 


with  the  equality  sign  only  when  ki  =  k2-  It  is  then  a  consequence  of  (3.28)  that  both  the 
real  and  the  imaginary  parts  of  the  zero  of  i?(^,  y)  in  the  ^  plane  are  positive  when  /ci  >  /c2, 
at  least  for  small  values  of  y. 


IV.  INTEGRAL  EQUATION  FOR  THE  CURRENT 

In  the  preceding  section,  the  function  E{^,y),  the  Fourier  transform  of  the  electric  field 
Ex  except  for  an  overall  constant  factor,  is  found  to  have  a  zero  at  .^  =  /cq  say,  where  ko  is 
given  approximately  by  the  square  root  (with  positive  real  part)  of  the  right-hand  side  of 
(3.28).  That  there  is  one  and  only  one  such  ko  plays  a  central  role  in  formulating  an  integral 
equation  for  the  current  on  the  linear  antenna  at  the  boundary  of  two  dielectric  materials. 
Of  course,  ^  =  —ko  is  also  a  zero  of  E{^,y). 

This  ko  has  a  number  of  important  properties.  First,  as  already  mentioned,  in  the  limiting 
case  where  ki  =  /c2,  this  ko  is  equal  to  ki.  The  details  of  how  this  limit  is  approached  is 
to  be  found  in  Appendix  C.  In  this  limiting  case,  this  ko  =  ki  is  essential  in  obtaining  the 
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Im  ^ 


-k  -k 


0 


k- 


^-plane 


Re^ 


FIG.  5:  A  usual  choice  of  the  hrst  sheet  in  the  ^-plane.  With  this  and  a  number  of  similar 
choices,  the  zeros  of  y)  at  ^  =  ko  and  ^  =  —ko  are  on  the  hrst  sheet. 


Hallen  integral  equation.  Secondly,  in  the  present,  more  general  case  of  two  media,  this 
ko  can  be  used  in  an  analogous  manner  to  give  an  integral  equation  for  the  current  on  the 
linear  antenna.  This  is  to  be  discussed  in  the  remainder  of  this  section.  The  importance  of 
this  ko  is  the  underlying  reason  why  it  is  studied  in  such  detail  in  Sec.  III.  Thirdly,  as  seen 
from  (3.30),  this  ko  has  a  positive  imaginary  part.  That  ko  is  not  real  even  in  the  absence 
of  dissipation  leads  to  features  that  are  not  seen  in  the  usual  case  of  the  linear  antenna  in 
a  uniform  medium.  Fourthly,  this  ko  is  in  the  hrst  sheet  of  the  complex  ^-plane  as  usually 
chosen — see  Fig.  5.  If  this  zero  of  E{^,y)  had  been  on  the  second  sheet  in  a  way  similar  to 
the  Sommerfeld  pole,  the  linear  antenna  on  an  interface  would  have  qualitatively  diherent 
properties.  Dehne 

=  (4-1) 

Then,  similar  to  E{^,y),  G  is  a  function  of  ki,  k2,  ^  and  y.  It  follows  from  (3.4)  and  (3.1) 


that 


(4.2) 


where 


/OO  poo  jp  poo 

e*«^G'(^,  y)=  ^2  _  1.2  /  e 

■OO  J  —00  S  ^0  J  — OO 


Mx+r,y)  7l(fe|-^^)+72(fe^-e^) 


MN 


(4.3) 
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is  a  function  of  /ci,  /c2,  x  and  y. 

In  order  to  gain  some  physical  intuition  for  this  G,  consider  once  more  the  special  case 
of  ki  =  k2,  where  this  G  simplifies  to 


-I  /»oo  /»oo  -I 


(4.4) 


This  integral  can  be  easily  evaluated  using  polar  coordinates  to  give 


G{x,y) 


-7ri 

\/x‘^  +  y^ 


(4.5) 


Therefore,  in  this  special  case,  G  can  be  used  to  operate  on  the  current  to  give  the  vector 
potential.  In  the  present  case  of  two  media,  although  no  vector  potential  can  be  usefully 
defined,  this  G  of  Eq.  (4.4)  plays  the  role  of  going  from  the  current  on  the  antenna  to  a 
“vector  potential.”  Thus  G  can  be  used  to  give  an  integral  equation  for  the  current. 

At  the  beginning  of  Sec.  II,  it  is  noted  that  a  thin  strip  antenna  of  width  4a  is  electro- 
magnetically  essentially  the  same  as  a  circular  one  of  radius  a.  Since  the  usual  argument  for 
this  equivalence  is  not  modified  for  the  present,  more  general  situation  of  two  media,  it  will 
not  be  repeated  here.  This  leads  to  the  integro-differential  equation 

+  [  dx' K{x  -  x')I{x)  = —E^p{x)  (4.6) 

\dx^  J  cn/io 

as  a  consequence  of  (4.2)  and  (4.3).  Here  the  kernel  K  of  the  equation  is 


K{x)  =  —  [  dO  G{ki,k2,x,2asm^6),  (4.7) 

271-  J-n 

I  is  the  current  on  the  antenna,  and  E^p{x)  E^^^{x,  0)  is  the  x-component  of  the  external 
electric  field  applied  to  the  thin  strip,  evaluated  at  y  =  0. 

Figure  6,  obtained  from  the  numerical  work  of  Sec.  VI  and  the  appendices,  shows  the 
kernel  K{x)  for  the  case  ki/k2  =  2  and  k2h  =  5,  for  two  values  of  a/h  =  10“^  and  a/h  = 
10“®.  Unlike  the  wire  in  a  homogeneous  medium,  the  kernels  for  the  two  cases  of  a/h  differ 
noticeably  even  for  values  of  a;  ^  a. 

It  is  also  of  interest  to  compare  the  kernel  K{x)  for  the  case  above  with  a/h  =  10“^  to 
the  Hallen  kernel  in  [1]  for  a  single  homogeneous  medium  adjusted  to  have  a  propagation 
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constant  /cq,  and  scaled  to  have  the  same  limiting  behavior  as  K{x)  as  x  tends  to  zero.  This 
Hallen  kernel  is 


def  1  exp(iA:o\/a;2  +  ‘^o?  sin^(6'/2)) 

{kihy  +  {k2hy  Jq  -s/ x"^  +  sin^ (^/2) 

Figure  7  compares  K{x)  against  this  H{x)  for  the  case  ki/k2  =  2  and  k2h 

a/h  =  10“'^. 


(4.8) 


5,  with 


V.  METHOD  OF  DETERMINING  RESONANCES 


We  believe  that  the  resonant  freqnencies  for  scattering  by  a  wire  are  relatively  independent 
of  the  incident  held,  regardless  of  whether  the  wire  in  a  homogeneons  medinm  or,  as  is  the 
case  here,  at  an  interface.  If  there  is  no  incident  held,  that  is,  if  Eg_p{x)  =  0  on  —h  <  x  <  h, 
then  Eq.  (4.6)  specializes  to 

J  =  0-  (5-1) 

At  hrst  glance,  one  might  expect  the  only  solntion  to  be  I{x)  =  0.  However,  the  key  to 
dehning  resonance  is  to  note  that  the  kernel  K  depends  not  only  on  position  along  the 


FIG.  6:  Two  cases  of  the  kernel  when  ki/k2  =  2  and  k2h  =  5;  a//i  =  10  ^  and  a//i  =  10 
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FIG.  7:  Kernel  K{x)  compared  with  Hallen  kernel  H{x)  for  case  of  k2h  =  5,  ki/k2  =  2,  and 
a/h  =  10““^. 

wire,  but  also  on  the  angular  frequency  a;,  on  which  ki  and  k2  depend  through  Eq.  (3.2). 
Therefore  the  solution  to  the  integral  equation  (5.1)  depends  on  uj  and  for  certain  discrete 
special  values  of  cu,  cui,  uj2,  . . . ,  the  equation  has  non-zero  solutions.  These  values  are 
expected  to  be  complex,  so  that,  listed  in  increasing  order,  we  have 

RecUn  =  the  n-th  resonance  frequency, 

— Ima;^  =  the  width  of  the  n-th  resonance.  (5.2) 

The  real  parts  of  the  complex-valued  solutions  are  the  resonant  frequencies,  while  the  imag¬ 
inary  parts  give  the  half-width  at  half  height  for  each  resonance,  under  the  assumption  that 
the  resonant  behavior  for  the  n-th  resonance  is  proportional  to  l/(a;  —  a;„).  The  task  is  to 
hnd  the  approximate  values  of  the  complex  ujn- 

To  determine  the  special  frequencies,  we  use  a  Galerkin  method,  and  for  this  we  need  to 
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symmetrize  Eq.  (4.6).  For  that,  noting  the  boundary  condition  /(±h)  =  0,  we  observe  that 


d  r  r 

—  /  dx'  K{x  —  x)I{x')  =  /  dx' 


dx 


'-h 


J  —h 

=  [  dx' 

J-h 

rh 


'-h 


^K(x-x')  I(x') 
OX  J 

-^K{x-x')  I{x') 

d 


=  I  dx' ^ —-^[K{x  —  x')I{x')]  +  K{x  —  x') 


=  -K{x-x')Hx%z''_,  +  j  dx'K(x-x')^,I(x') 

d 

=  /  dx  Kix  —  x')——I{x'). 

J-h  ox' 

With  the  use  of  this  equation,  Eq.  (5.1)  takes  the  symmetrical  form 


d 


(5.3) 


[  dx' K(x  —  x')-^I(x')  +  kl  [  dx' K(x  —  x')I(x')  =  0. 
ux  J —fi  ux 


(5.4) 


’-h 


This  is  the  form  of  the  Pocklington  integral  equation  to  which  we  shall  apply  Galerkin’s 
method.  We  deal  with  the  resonances  that  have  symmetric  currents  separately  from  reso¬ 
nances  for  which  the  current  is  antisymmetric  in  x.  For  a  trial  current  /Gai(a()  in  the  Galerkin 
subspace,  after  an  integration  by  parts  Eq.  (5.4)  takes  the  approximate  form: 


d 

/  dx 
J-h 

I  dx' K{x  -  x')-^Igixi{x') 

rh 


+  kl  dxlcaiix)  /  dx' K{x  -  x')Ig^i{x')  =  0- 


(5.5) 


’-h 


’-h 


1.  Galerkin  subspace  for  resonances  with  symmetric  and  antisymmetric  currents 


For  sufficiently  small  a,  one  expects  resonant  frequencies  at  complex  values  of  ko  near 


def 


(5.6) 


for  n  =  1,  2, . . . ,  with  the  currents  corresponding  to  odd  values  of  n  symmetric  in  x  while  the 
currents  corresponding  to  even  values  of  n  are  antisymmetric  in  x.  In  applying  a  Galerkin 
method,  a  critical  question  is  this:  should  the  Galerkin  subspace  depend  on  the  frequency  uj 
or  not?  For  the  symmetric  case,  there  are  good  reasons  to  believe  that  the  resonance  current 
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is  given  roughly  by  the  shifted  cosine  form  coskx  —  cos  kh  [1].  One  might  expect  to  get  a 
good  approximation  using  a  Galerkin  subspace  of  dimension  1  by  choosing  the  current  to  be 
the  shifted-cosine;  however,  use  of  the  shifted-cosine  results  in  possibly  spurious  values  of 
the  resonant  frequency.  Examination  of  how  these  values  arise  shows  that  they  are  indeed 
spurious  artifacts  of  the  shifted-cosine  form;  hence  we  need  to  attend  more  carefully  to  the 
choice  of  the  Galerkin  subspace.  The  shifted-cosine  form  leads  to  a  Galerkin  subspace  de¬ 
pendent  on  the  frequency  cu  ~  A;2-  Question:  should  the  Galerkin  subspace  depend  on  the 
frequency  u  or  not?  Having  encountered  spurious  values  from  the  shifted-cosine  form,  we 
choose  instead  a  Galerkin  subspace  that  is  independent  of  the  frequency  oo.  This  indepen¬ 
dence  essentially  determines  the  Galerkin  subspace  of  dimension  1.  To  begin  with,  it  can 
depend  only  on  the  geometric  parameters  h  and  a.  Although  we  could  force  the  Galerkin 
subspace  to  depend  on  a,  we  cannot  see  how  to  do  this  in  a  physically  sensible  way.  Thus 
we  take  the  Galerkin  subspace  to  depend  only  on  h.  Then  we  can  hardly  avoid  choosing 
~  cos  K.nX  =  cosinnx /2h)  for  the  symmetric  case  where  n  is  odd.  Gorrespond- 
ingly,  for  the  cases  of  resonance  in  which  the  current  is  antisymmetric  about  the  center 
point  of  the  wire,  for  which  n  is  even,  our  Galerkin  subspace  consists  of  the  single  function 


rasym 

-'Gal 


sin  KnX. 


2.  Resonances 
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and  carry  out  the  differentiations  in  Eq.  (5.7)  to  obtain  for  n  odd: 


-K-n^sn  +  k^Icn  —  0, 


(5.10) 


which  in  a  form  more  convenient  for  calculation  becomes, 


~  kQ){Isn  +  Icn)  +  (^n  +  kQ){Isn  ~  -^cn)  —  0; 


(5.11) 


Similarly  for  the  resonances  with  antisymmetric  currents  so  that  n  is  even  one  finds: 


ch 

/  dx 
J-h 

■  d 

Sm  KnX 

Ox 

fh 

+k^ 

/  dx  sin  KnX 

J-h 

J  —h 

[  d 

'-h 


^  /  /  <9 

dx'  Kix  —  x)  — — -  sin  Knx' 
ox' 


(5.12) 


which  implies 


(^n  ~  k‘^){Isn  +  Icn)  ~  (^n  +  k‘^){Isn  ~  Icn)  —  0, 


(5.13) 


We  now  make  a  zero-th  order  check.  On  replacing  K{x  —  x')  by  a  delta  function  5{x  —  x') 
and  carrying  out  the  integrals  in  Eq.  (5.7)  and  using  the  definition  of  Kn  =  {n7i/2h),  the 
result  for  n  odd  is  +  /cp  =  0;  carrying  out  the  same  procedure  on  Eq.  (5.12)  yields 
this  relation  for  n  even;  confirming  our  claim  that  for  sufficiently  small  radius,  the  resonant 
propagation  constants  should  be  near  K,n- 

Returning  to  deal  with  K{x  —  x')  and  not  just  the  delta  function,  one  holds  fixed  the 
geometrical  and  material  parameters  h,  a,  ei,  62  and  /io  while  varying  the  angular  frequency 
uj  in  order  to  find  complex  values  of  u  for  which  Eq.  (5.11)  for  the  symmetric  resonances 
and  Eq.  (5.13)  for  the  antisymmetric  resonances  have  non-zero  solutions.  Varying  uj  in  the 
complex  plane  implies  also  varying  the  propagation  constants  /ci,  /c2  via  kj  =  for 

j  =  1,2;  and  the  dependence  of  ki  and  k2  on  (complex- valued)  oj  varies  with  ko  as  determined 
by  Eq.  (3.28). 

The  next  step,  derived  in  Appendix  D,  is  to  reduce  the  double  integrals  to  single  integrals 
to  obtain,  for  all  resonances,  regardless  of  whether  n  is  even  or  odd. 


i^n  kQ)lY,n  —  (^n  +  ^o)9An; 


(5.14) 
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where  we  define 


^sn  “t“  ^cn  ^  /  dxl^{x^{^h  X^COSK^X^ 

Jo 

,  f  2 

Iau  =  (-l)”(/sn  -  /cn)  =  -  /  dxK{x)smKnX. 

Jo 


(5.15) 


3.  Scaling  to  h  =  1 


For  calculational  purposes,  we  save  effort  by  recognizing  that  the  resonant  frequencies 
scale  with  the  wire  half-length  h.  Thus  one  can  set  h  to  1  and  obtain  resonant  frequencies 
for  other  values  of  h  by  dividing  by  h.  With  h  =  1  we  have  =  mi/2  and  the  task  is  to 
solve  the  equation  for  resonances,  namely,  regardless  of  whether  n  is  even  or  odd. 


\m\l2f  -  kl]I^n  =  [(W2)^  +  kl]Ij^n, 


(5.16) 


where  with  h  =  1  and  Kn  =  n'x/2  we  have 


-fsn  *=  Isn  +  Icn  =  ‘J  /  dx  K  (x)  {2  -  x)  COS  , 

Jo  ^ 

Iau  =  (-1)  {Isn-Icn)  =  —  /  dx  K  (x)  sm  . 

mr  /n  2 


(5.17) 


Similarly  to  the  case  of  the  wire  in  a  homogeneous  medium,  when  the  strip  width  4a  is  much 
smaller  than  all  the  other  dimensions  in  this  scattering  problem,  we  can  approximate  the 
kernel  K{x)  (defined  in  Eq.  (4.7))  by  G{x,a)  defined  in  Eq.  (4.3),  so  that  we  have 


K{x)  ^  G{ki,k2,x,a)  =  /  -r - /  drie 

J-oo  e  J-oo 


ii^x+va)  7l(fe|  +72(^1  -C^) 

MN  ■  ^ 


We  will  use  this  approximation  to  determine  resonances. 


VI.  NUMERICAL  ANALYSIS  OF  RESONANCES 


Here  we  address  the  numerical  analysis  of  Eqs.  (5.16)  through  (5.18)  for  a  wire  scaled 
to  a  half-length  of  h  =  1.  Recall  that  ko  is  the  zero  of  the  Fourier  transform  of  the  kernel 
K{  x).  Later  we  will  use  the  abbreviation 

/c4  =  kik2 !  \Jkl  +  kl  =  k2J  (6-1) 
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so  that  for  the  relative  dielectric  constant  e  real  and  positive,  ki,  /c2,  and  /C4  as  complex 
parameters  all  have  the  same  argument. 


For  keeping  track  of  branch  cuts,  in  this  section  we  study  imaginary  values  of  these 
parameters,  and  most  of  the  following  will  deal  in  “kappa’s”  and  corresponding  7’s  defined 


by 

(For  j  =  0, 1,  2, 4)  kj  =  inj, 

(6.2) 

(Forj  =  l,2) 

(6.3) 

where 

7j  +  7^  +  1^]- 

(6.4) 

For  M  and  N  defined 

in  Eqs.  (2.5)  we  have 

M  =  i(7i  72),  and  N  =  +  ^172). 

(6.5) 

[Note:  k’s  as  defined  here  have  nothing  to  do  with  use  of  the  Kn  symbol  in  the  sections 
above.]  The  advantage  to  this  notation  is  that  for  real  k’s,  the  branch  points  and  branch 
lines  present  no  trouble.  We  then  determine  the  proper  branches  involved  in  various  integrals 
by  following  a  path  from  positive  real  k’s  to  whatever  complex  values  we  want  with  the  rule 
of  “no  cross  branch  lines”  which  tells  us  how  to  deform  branch  lines  so  as  to  avoid  such 


crossings.  Changing  from  k's  to  k's  and  recognizing  that  numerical  integration  over  the 


kernel  K{x)  takes  almost  a  minute,  we  put  Eq.  (5.16)  into  the  form 


/h/i  \“  9  /  \  luiju  ^  /  iiJi  \“  2 

—  +Kq  {2-x)  cos  — - —  —  -Ko  sm 


K{x)  =  0.  (6.6) 


To  express  the  kernel  K{x)  given  in  Eq.  (5.18),  we  define 

Tf  =^7i  +  72  and  N  =  +  «i72- 


(6.7) 


After  dropping  constants  that  cancel  out  in  the  equations  above  for  resonant  frequencies. 


for  real-valued  frequencies  the  kernel  of  the  integral  equation  can  be  taken  to  be: 

.i(^x+na)  7l('C^  +  +  72('C^  +  K-l) 


K  x)  = 


/_oo  7-00 


fin  (=7?^+^“)  7i(C  +  ^2)  +72(C  +  K-i) 
(71 +  72)(«i7i  +  «?72) 


1-00  +  ^0  7-00 


^i(^x+rja) 


1 

^ - 1 - ^ 

M  N 
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Remark:  for  complex-valued  resonant  frequencies  of  interest,  however,  the  in¬ 
tegration  path  must  be  deformed  to  avoid  crossing  branch  cuts.  As  a  result  the 
expressions  to  be  developed  for  the  kernel  will  be  formulated  in  two  steps:  (1)  in  accordance 
with  Eq.  (6.9),  and  (2)  by  analytic  continuation  for  the  range  of  resonant  frequencies  that 
have  positive  real  parts  but  negative  imaginary  parts,  corresponding  to  K2  =  where 

a  exceeds  n/2. 

The  kernel  K{x)  depends  on  the  resonant  frequency  through  its  dependence  on  ki,  K2, 
and  Kq,  which  itself  is  a  function  of  Ki  and  K2,  provided  by  the  MATLAB  program  f  indkOA.m 
in  Appendix  F. 

The  integrand  has  poles  at  .^  =  ±mo.  Adapting  Eqs.  (3.28)  and  (3.29)  from  ko  to  kq, 
one  shows  that  for  a  real-valued  dielectric  constant  e  =  4,  over  the  range  0  <  a  <  1.8, 
has  a  relatively  large  positive  real  part  and  a  small  positive  imaginary  part.  It  is  this  pole 
that  contributes  in  later  work  when  we  deform  the  integration  path. 


1.  Summary  of  the  kernel  transformed  for  numerical  calculation 

The  double  integrals  that  define  the  kernel  need  to  be  put  into  a  form  that  appears  much 
more  complicated  but  is  more  amenable  to  numerical  work.  Define 


J^{u)  ^  [(1  K2u)e  -  (1  +  K,iu)e  . 

o  /O  L 


u 


(6.10) 


Then  we  find  below  that 


K(x)  = 


Att 


Kl  —  Ko 


T[^/x^  +  a2) 


0^2  _  ^2  _  2  j-x  _ 

H - ^ ^  /  dx  sinhKo(a;  —  x')Ti(\/x''^  -t-  a^) 

+  K2{Q)  cosh  KqX  +  K:ii{x)  +  Ks2ix), 


(6.11) 


where  as  shown  in  Appendix  E, 


K2{0)  = 


TT 


—  Ko 


2kI 


k\  +  kI 
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+ 


K^a 
Ki  —  K2  -\- 
TT 


ni  -  Kg  Ki  K2  -  “  y  (  « V 


Kl  —  Ko 


2ko 
2^0 

kI  +  k| 


In 


K2  -  Ko 


K2  -  Kl 


2Kn 


In 


Kl  -  Ko 


(6.12) 


KT  —  Ko 


Ko 


('ni+7-i) 


ff2  _  ^2  k"^  - 

+  ^ ^  ln(K2  +  Ko)  -  ^ ^  ln(Ki  +  Ko) 


+  0{a^  Ina). 


(6.13) 


Ko  Ko 

In  the  equation  above,  Ki  is  a  modified  Bessel  function  [4].  For  the  remaining  parts  of  the 
kernel  we  have 


K^iix)  = 


n 


Kl  —  Ko  \  X 


k2  _  k2  +  ^  +  ^  ) 

X 


2  2  I  2K2  •  . 

Kg  -  Kl  + - +  ^  e 

X  x^  ' 


^  —  K,2X 


Ko(ki  +  kI-  Kg)  (^^T3i2(a;)j 


+ 


TTKq 

Kl  +  K2 


(6.14) 


with 

Ko, 


—Tsuix)  =  e 

TT 


In 


_  ^2 
Ko  K2 


+  iTT  +  e 


□  -Koa: 


(F'i[(ki  -  Ko)a;]  -  £'i[(k2  -  Ko)a;]) 


Ko2{x)  = 


—  e'^°^  (£'i[(ki  +  Ko)a;]  —  £'i[(k2  +  Ko)a;]) ; 

iKiK2K4 


4  4 

~  ^2 


(1^2(2^)  ~  1^1(2^))) 


(6.15) 


where  for  j  =  1,  2  we  define  Rj  by  Ri  ki/k2  and  R2  *=  K2/K1,  and 


Gj{x)  =  inKoe 


—  2Tii 


—  KQX  _ 


:  In 


Kq  —  k\  —  iK^Rj 


\/kI-  kI  a/k|  -  k^  +  iK^Rj 

"V^-i  ((  +  i)2exp[-K4(C  +  l)a:] 


dC 


(6.16) 


•^0  [(c  + 1)^  —  ki/ K4]  \/  c(c  + 1) 

Note  that  the  correct  branch  of  the  exponential  integral  £'i[(k2  —  Ko)a;]  in  the  expression 
for  K3i(a;)  requires  adding  27ii  in  case  the  imaginary  part  of  K2  —  Kq  is  positive.  Other 
important  details  concerning  the  branches  of  the  above  functions  are  given  below.  For 
X  <  0.001  we  make  use  of  the  relations 
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2 

X 


9  9  2/x'l  2 

«0  -  «2  +  -  +  ^ 


^  4  1 


X 


x^ 


=  2{ki  -  K2) 


2{k\  +  kI)  -  K1K2 


-  Kn 


+  -  «^2) 


^2  ,  21 
2  l%i  ^  l%2 
^0 


x  +  0{x^),  (6.17) 


Ei{ax)  —  Eiipx)  =  log(6)  —  log(a)  +  0{x‘^).  (6.18) 

Checks  made  on  the  correctness  of  these  equations  include  the  following:  We  checked 
that  for  a  =  7r/2  the  kernel  K{x)  as  computed  here  decreases  toward  zero  for  large  values 
x;  this  rules  out  any  large  error  in  7^2(0).  We  checked  that  for  a  relative  dielectric  constant 
e  near  1,  the  computed  kernel  approaches  the  free-space  case.  We  checked  the  K{x)  varies 
smoothly,  without  jumps,  as  kq  is  artihcially  scaled  by  a  parameter  that  ranges  from  0.3  to 
1.1,  and  we  checked  that  K{x)  varies  smoothly  as  a  ranges  from  0  to  2. 

2.  Putting  the  kernel  in  a  form  suitable  for  numerical  analysis 

We  rearrange  the  integrand  of  Eq.  (6.9)  into  three  pieces  in  a  way  that  makes  it  easy  to 
calculate  the  piece  most  sharply  peaked  at  a;  =  0 


^  ^  ^ 


2  2 

0  ^ 


M  '  N  + 


+  1- 


2«:2 


Ef  +  kI 


e 


1 

M  kI  +  kI 


kI  +  ^2 
N 


2 

M 


Correspondingly  we  have 


(6.19) 


K{x)  =  Ki{x)  +  K2{x)  +  Ksix) 


(6.20) 


where 


Ki{x)  = 
K2{x)  = 
K^ix)  = 


kI  +  kI 

2kI 


dC 


Ef  +  K-l 


’  —00 

POO 


^  m’ 
dC 


^  m’ 
i(^x+va)  ( ^4  +  4  _  A 


1-00  4  +  4  J-C 


4  +  4  J- 


dc 


e 


dT]E 


4  +  4  J-c 


e  +  4J-oo  '  \  N  M 

ac  r  ^  f4+4 

C2  ,  2  / 

r  +  4  J-00 


2 

M 


(6.21) 

(6.22) 


(6.23) 


''0  a  —00  \  N 

where  the  last  relation  follows  because  the  integrand  falls  off  fast  enough  as  r]  becomes 
large  to  make  K^{x)  smooth  as  a;  0,  so  that  in  the  integrand  we  replace  a  by  0  without 
signihcant  loss  of  accuracy.  [  Rechecked  and  found  correct  18  July  2009.  ] 
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APPENDIX  A:  COMPARISON  WITH  HOMOGENEOUS  MEDIUM 


In  contrast  to  the  assumption  in  the  main  report  that  ki  >  k2,  in  this  appendix  we  discuss 
the  case  ki  =  k2  =  k,  for  which  Eq.  (2.4)  becomes 


Ex{x,y)  = 


cn/io 

47J-2 


where 


E,{C,y)={k^-e)  / 


(Al) 


(A2) 


2k‘^sjk‘^  — 

is  the  Fourier  transform  from  a;  to  of  Ex{x,  y).  From  the  factor  (/c^  —  in  Ex{C,  y)  follows 
the  relation  for  Ex{x,y): 


EA^.v)  =  + 


47r2  \dx'^ 

{  d‘^  ,  ,2 


'  — OO 

noo  poo 


(u/io  ^2 1  / 


'  — OO  j  — OO 


2k‘^\/k‘^  — 


47r2  \dx'^ 


— +  /c  Go{x,y), 


where 


G„{x.y)=  e-'-v^. 


(A3) 


(A4) 


This  Go  generates  the  kernel  of  the  integral  equation  of  the  Pocklington  type  for  the  current 
for  the  scattering  by  a  thin  perfectly  conducting  strip: 


poo  p2a 


E"^^{x,y)  =  ^^(^  +  k‘^\  j  dx'  I  dy' G{x  -  x\y  -  y)J^{x\y'),  (A5) 

\  /  t/  — OO  t/  — 2r 

where  by  virtue  of  the  slimness  of  the  strip  {a  h  and  ka  1),  we  ignore  Jy,  the  current 
density  transverse  to  the  strip  axis  [1]. 
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APPENDIX  B:  FORM  OF  SERIES  EXPANSION  FOR  Ei(^,2/) 


In  this  appendix,  a  proof  is  given  that,  for  small  y,  the  series  expansion  for  the  Fi{y)  = 
Fi{^,y)  dehned  by  Eq.  (3.4)  is  of  the  form  of  (3.9),  i.e.. 


Fi  {y)  =  X]  +  «n]  • 


(Bl) 


n=0 


In  view  of  the  appearance  of  the  logarithm,  a  Mellin  transform  should  be  used.  Dehne 


n(C)  =  /  dyy-'*<Ffy). 


(B2) 


By  Eq.  (3.4),  this  is 


Fi(C)  =  2  dy 


Since 


Eq.  (B3)  reduces  to 


dyy 


0  ”  '  kly2  +  k^7i  Jo 

cosyy  = 


dyy 


-I+C 


cos  yy. 


r(i) 

r(¥)’ 


(B3) 


(B4) 


r  fi) 


^i72  +  ^271 


(B5) 


From  the  factor  r((C/2),  this  Ei(C)  has  poles  at  negative  even  integer  values  (and  zero)  of 
(.  Furthermore,  for  large  y,  l/{kfy2  +  ^ibi)  can  be  expanded  into  a  series  in  odd  powers  of 
1/y.  Using  analytic  continuation  as  usual  in  applying  the  method  of  the  Mellin  transform, 
the  structure  of  this  series  leads  to  an  additional  factor  of  1/  (C  +  2n)  near  negative  even 
integer  values  of  (,  but  no  such  corresponding  factor  for  negative  odd  integer  values  of  (. 
This  means  that  Ei((C)  has  double  poles  at  negative  even  integer  values  of  but  is  analytic 
at  negative  odd  values  of  (.  Therefore 


—  ttr 


+ 


(B6) 


;^L(C  +  2n)2  <;  +  2nJ 

gives  explicitly  the  positions  and  nature  of  the  poles  in  the  left  half  of  the  C-plane-  The 
desired  result  Eq.  (Bl)  follows  from  Eq.  (B6)  by  inverting  the  Mellin  transform  (B2). 


Note  that  we  have  paid  no  attention  to  the  issue  of  whether  the  series  on  the  right-hand 
side  of  Eq.  (Bl)  is  convergent  or  not.  For  the  present  purpose  of  studying  the  properties 
of  a  linear  antenna  on  an  interface,  only  the  terms  oq  and  Oq  are  of  interest,  and  thus  the 
convergence  of  the  series  is  irrelevant. 


APPENDIX  C:  CASE  OF  fci  CLOSE  TO  fca 


In  this  appendix,  the  case  in  which  ki  and  /c2  are  close  to  each  other  is  discussed. 

When  /ci  =  /c2,  the  zero  of  E  is  given  hy  ^  =  ki  (=  /C2).  This  limiting  value  is  exact  for  all 
values  of  y.  Nevertheless,  only  the  case  of  small  y  as  defined  by  (3.7)  is  to  be  treated.  This 
special  case  corresponds  to  the  wire  on  the  interface  being  thin,  and  is  especially  simple  to 
deal  with,  because  the  explicit  results  (3.28)  and  (3.29)  of  Sec.  Ill  can  be  used. 

Since  ki  and  /c2  are  close  to  each  other,  it  is  convenient  to  define  k  >  0  and  e  >  0  so  that 


kl  =  k'^  +  e 


^  and  kl  =  k‘^ 


(Cl) 


With  e  /c,  various  quantities  are  expanded  in  the  small  parameter  jk^ . 
Consider  first  the  real  part  on  the  right-hand  side  of  Eq.  (3.29): 


TT  1 

U2  1  u2  u2u2 

TT  1 

r  _  4  1 

37r 

0  U2  _  U2 

Zj  rv‘2 

2  y{kt  +  kt_)/2 

^  2  ^ 

Vk‘^  +  e'^_ 

(C2) 

t-  L  v/c"*  +  6“*]  O  K- 

This  is  indeed  positive,  consistent  with  Eq.  (3.30). 

The  imaginary  part  on  the  right-hand  side  of  Eq.  (3.29)  can  be  treated  in  a  similar  way, 
only  slightly  more  complicated: 


1 

'X  + 


k\^/2  -\-  \/kf  -[-  /cl 
~  ^2  \/ {kf  k^)/^  -f 


klk. 


2U2 
2 


In 


In 


k"^  + 


_  _1  k^-€^  1 

2  26^  ypT|r^  ^2  _  g2  _|_  _|_  g4 


2  ^  2£2  V  2  P 


1 


29 


rsj 


rsj 


rsj 


e^\  2PJ  +  3\2PJ 

^  _  1  ^ 

2P  8p)  ^  24  k\ 


7  e^ 

8  ¥' 


(C3) 


In  particular,  this  calculation  shows  the  important  role  played  by  the  hrst  term,  — 

Finally,  the  substitution  of  Eqs.  (C2)  and  (C3)  into  (3.28)  gives  the  real  and  imaginary 
parts  of  p,  where  ^  is  the  zero  of  y): 


~  P 


7  1 

8  P  ln(|/e) 


(C4) 


and 


fm4  ~ - 


8  HyTY 

where  the  terms  7  — ln2  — iTr/d  in  the  denominator  of  Eq.  (3.28)  have  been  neglected  because 
they  are  much  smaller  than  |  ln(|/e)|. 

The  results  (C4)  and  (C5)  are  believed  to  be  valid  without  the  condition  py  1  of  Eq. 
(3.7). 


APPENDIX  D;  INTEGRATIONS  EOR  !<.„  AND  J,„ 


Because  K{x  —  x')  depends  only  on  the  difference  between  x  and  x',  it  must  be  possible 
to  rewrite  the  integrals 


as  single  integrals, 
argument.  Dehne 


Isn  /  dxsintYnX  I  dx'K{x  —  x')  sin  KnX, 

J-h  J-h 

/h  i*h 

dx  cos  KnX  /  dx'K{x  —  x')  cos  KnX 

■h  J-h 


f*h 


(Dl) 

(D2) 


We  do  this  without  using  the  fact  that  K  is  an  even  function  of  its 


/h  ph 

dx  /  dx'  K{x  —  x)e 

-h  J-h 


/\  Aol  X 


(D3) 
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Changing  integration  variables  x  — ^  —x'  and  x'  — ^  —x  produces  the  relation 


Ie{  —  Ct',  —a)  =  Ie{ct,  a'). 

(D4) 

With  this  relation  one  expresses  Isn  and  Icn  as 

Isn 

^  [2-/^e(^n:  ^n)  ^n)  -^e( 

(D5) 

Icn 

“  Kfjj  +  l^n)  “1“  -^e(  ^n)]: 

(DO) 

SO  that  we  have 

Isn  “1“  Icn  2  ^n)  “1“  -^e(  ^n-j  ^n)]: 

(D7) 

Isn  Icn  lei^^n')  ^n)* 

(D8) 

For  the  reduction  to  single  integrals  we  compute 


Je  («,«') 


/h  ph 

dx'  / 

-h  J-h 

/h  ph—x' 

dx'  / 

■h  J —h—x' 

/o  ph  p2h  r^—y 

dyK{y)e^^y  /  +  /  dy  K{y)e^^y  / 

■2h  J —h—y  Jo  J —h 

t{a  +  a')lJ_2h  V  J 

+  f  dy  K(y)e^^y  ( eh“+«')(^-y)  _  ^-i{»+a')h'\ 


+  J  dy  K{y)e^^y  j 

=  — ^ - -  [  [  dy  Kiy)  ( eh«+«')/igiay  _  ^-iia+a')h^-Wy'\ 

i{a  +  a')lJ_2h  V  / 

f2h  1 

The  case  «'  =  — o;  is  worked  out  directly  to  show 

/O  p2h 

dyK{y){2h  +  y)e^^y+  /  dy  K{y){2h  -  y)e^<^y . 
-2h  Jo 

Substitution  of  Eqs.  (D9)  and  (DIO)  into  Eqs.  (D7)  and  (D8)  yields 


(DIO) 


h^  +  Icn=  /  %A’(y)(2/i  +  y)cosK„y+  /  dy K{y){2h  -  y)cosK„y,  (Dll) 

J-2h  Jo 

^  ~  1^0  i*2h 

Isn-Icn  = - /  dy  K{y)si\i[Kn{2h  + y)]+  /  dy  K{y)si\i[Kn{2h  -  y)]  .  (D12) 

\_J -2h  Jo 
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y 


In  the  special  case,  which  we  have  here,  in  which  Kiy)  =  K[—y)  and  =  {n7i/2h),  these 
simplify  to 

hn  +  hn  =  2^  dy  K{y){2h-y)  COS 

Isn-Icn  =  (-1)”^^  dyK{y)sin^^^,  (D13) 

which  are  the  single  integrals  that  we  wanted  to  obtain. 


APPENDIX  E:  EVALUATION  OF  K^,  K2,  AND 


The  integral  in  Eq.  (6.21)  that  dehnes  Ki  is  evaluated  exactly  as  follows: 


o  /»00  POO  -I 

Ki{x)  =  .  /  di  / 

kI  +  kIJ_^  7-00  M 


^00  POO 


/  drje 


i[^x-\-r}a) 


"1“  ^2  J  —00  J  —00 


\/rf  +  kl  +  \/rf  +  kl 


Let 


p  +  a^. 


and  change  to  polar  coordinates 


X  =  pcoscj),  a  =  psin0,  .^  =  Acosa,  and  r/  =  AsinQ;, 


(El) 

(E2) 


(E3) 


(E4) 
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leading  to 


.2  I  ,^2 


2  ^27r  ^oo  gipAcos(^— a) 

.,  ,  da  XdX — ,  - , 

“1“  ^2  Jo  Jo  +  A^ 


47r  /■“ 

2  /  AdA 


Jo(Ap) 


«:?  +  do 


\/ K, I  -\-  X‘^  +  Y^K2  +  A^ 


47r 


KT  —  Ko 


-[J(p,fi;i)  -X(p,  Ka)], 


(E5) 

(E6) 

(E7) 


where  we  define 


I{p,k)  [  XdXJo{Xp)  ^ 


=  \  K- - —^ - —  /  XdX  ,  = 

dp^  pdpj  Jq  ^/k‘^  +  A2 

(9^  1  (9  A  e“'^^ 


=  \  K  - 


(9p2  pdp)  p 


=  — ^(l  +  «p)e 


(E8) 

(E9) 

(ElO) 

(Ell) 


^2  2  (9 

in  which  we  used  the  relation  (  7—  H — - — h  A^  )  Jo(Xp)  =  0.  Putting  all  this  together  yields 

op^^  pop 


Ki(x)  = 


Ati 


(Kf  -  kDp^ 
where  x  is  related  to  p  by  Eq.  (E3). 

For  /cap  1,  it  follows  that 


[(1  + fi;ap)e  -  (1  + fi;ip)e 


(E12) 


TT  2  4  Kf  +  K1K2  +  1  2  .  2> 

Klix)  —  — ^ K - — - h  “(^1  +  ^2)p 

^  '  kI  +  kI[p  3  Ki  +  Ka  2^  ^  _ 


+  0{p^). 


(E13) 


1.  Evaluation  of  K^^x) 


From  Eqs.  (6.21)  and  (6.22)  it  follows  that 

-  X  j  K,(x)  =  k,(x).  (E14) 

Solving  the  differential  equation  and  imposing  the  symmetry  condition  that  K2{x)  = 
K2{—x)  yield 

(2  2  p 

Ko - /  dx' sinh[fi;o(a^  —  a;^)].f^i(a^0  + -^2(0)  cosh(fi;oa;).  (E15) 

2ko  J  Jo 
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The  coefficient  K2{0)  is  a  function  of  a  which  increases  without  bound  as  a  0;  it  is 
evaluated  as  follows  (and  implemented  in  K20.m).  By  the  dehnition  of  K2  in  Eq.  (6.22), 

/■“  dC  1 


1^2(0)  =  1 


1^1  +  1^2/  J-00  +  1^0  J-c 


drj  e 


irja 


Kt  —  Ko 


Kt  — 


2kI 


k\  +  K2 

k\  +  kI 


dr]  e 


irja 


7i  +  72 

dC 


drj  e 


iTja 


.71  72 


+ 


drj  e 


I'qa 


di 


Loo  L  +  L 


rf  +  k\-  kI  rf  +  kI-  kI 


7i 


72 


Kt  —  Ko 


2kI 


kI  +  K-l 


{T1+T2), 


where  we  dehne 


(E16) 


rji  del 

J 1  — 


drj  e 


trja 


.71  72 


a  ^ 


_  _27r(fi;i  -  K2), 


rji 

i  9  — 


(E17) 

(E18) 

(E19) 


where  for  j  =  1 ,  2  we  dehne 


Liv)  =  + 


dC 


LI  +  L) \jL  +  L  +  L] 


Kq 


=  —  JL  +  L  -  kI  In 


Ko  - 


«o  +  T  +  L  -  L 


Ko 


=  —  JL  +  L  -  kI  tt  +  i  In 


rf  -\-  L  -  Kq  +  iKo 


(E20) 


(E21) 


(E22) 


'rj^  +  Kj  -  Kq-  ikq 

The  last  form  is  convenient  for  obtaining  the  behavior  for  large  \r]\  [MATLAB  check  of 
Eq.  (E21)  against  Eq.  (E20)  passed.]  Thus  we  have 


T2  =  T21  +  T22, 


(E23) 


34 


where 


where 


+  0(a^  In  a). 


(E24) 

(E25) 

(E26) 


(E27) 

(E28) 


(In  the  above  equations,  unlike  the  body  of  the  paper,  Kq  and  Ki  denote  modified  Bessel 
functions  [4].) 

Equations  (E26)  and  (E28)  imply 


+  0(0^  In  a).  (E31) 


(MATLAB  checks  OK.) 
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For  T22  we  have 


„  tM  * 

J-22  —  — 

Kq 


dr]  e 


irja 


rf  +  k\  —  Kq  In 


^Jrf  +  kI-  kI  +  i 


IKq 


^Tj"^  +  kI-  kI-  mo 


r/2  +  —  Kq  In 


V^f  + 


^2  —  Kq  + 


\AF+ 


K: 


KO 


=  —  dr] 


rf  +  k\  —  Kq  In 


2-  k-q- 
^rj"^  +  kI-  kI  -  mo 


r/2  +  —  Kq  In 


^2  —  Kq  + 


\AF+ 


K: 


2-  t^O- 


=  TT  Ki  -  fi;2  + 


k; 


2  -  In  ^^2  +  «;o  K'i  -  f^o  +  f^o 


In 


(E32) 

(E33) 


2kq  K2  —  Ho  2fi;o  Hi  —  Ko_ 
where  the  last  expression  is  obtained  to  within  an  error  of  order  In  a  by  the  technique 
illustrated  in  Eq.  (E78).  That  expression  is  implemented  in  the  MATLAB  routine  T22.m. 
[May  need  to  modify  to  deal  with  branch.] 

Combining  Eqs.  (E31)  and  (E33)  yields 


To. 


n 


Hi  -  H2- 


^■2  _  ^2 
H2 

Ho 


\^1+1-^^  +  '^^H^2  +  Ho) 

Z  /  /t^Q 


K 


1  ^0 

Ho 


ln(fi;i  +  Ho) 


(E34) 


2.  Evaluation  of  K3 


From  Eq.  (6.23)  we  have 
Kofx)  ^ 


nj  +  hI 


die 


i^x 


+  ^0  J-00 


drj 


h\  +  h\ 
hI^I  +  h1^2 


7i  +  72 


(E35) 


This  can  be  rewritten  as 


Koix)  = 


hI  +  hI  J_, 


d(e'(^F3{0, 


(E36) 


where  we  define 


(  >=;  + 


Ho 


\hI'^i  +  hI'^2  71  +  72 


(E37) 
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drj 


2(71  -  72) 


[kI  -  +  Jo 

,  [2 ,  .2 ,  ^W2  \  f4i^i+e+v^)  ^i{4+e+v^)\ 

+  71  72+  J 


e 


(nl-nDi^l  +  e) 


f^i-4  ni  +  4+e,  nj+e 
2  2  °^«:i+e 


-\-G{Ki^  K,2,  0  ~  G{k,2,  ^1,  0 


In  (E39)  we  defined  G  by 

def  2  4 


drj 


'0  [{kI  +  ^i)  (7^  +  +  k^kI]  ’ 


_  i  K1K2K4,  -  in4i^2/ni 

2  \/^^  +  «!  \/C^  +  «!  +  iK4K2/Ki  ’ 

where,  as  follows  from  Eq.  (6.1), 

K1K2 


^4"^=  K2\/e/{e  +  1)  = 


y/Kl  +  K: 


Putting  this  together  yields 


where  for  0  <  a  <  n/2, 


Ks{x)  =  Ksiix)  +  Ks2{x), 


K34{x)  = 


1 


K^2ix)  =  i 


00 

£2  I  ^2 

2  — cxD  S  '  ^0 

00  pi^xf:2 

dC 


Kl  —  K. 


.  K1K2H4 


2  r  ^2  I  y2 

(^1  +  ^2  +  C  )  In 


2  -  («?  -  «^2) 


4  4 

~  ^2J-oo 


+  ^0  \/C^  +  K4 

^  I  \/C^  +  K\-iK4K2lKx  _  \/i^  iK4Ki/K2\ 

\/jJ  +  l^\  +  iK4K2/Ki  \/C^  +  «!  +  iK4Ki/K2  ) 


,  (E38) 


(E39) 


(E40) 

(E41) 

(E42) 

(E43) 

(E44) 


(E45) 


3.  Evaluation  of  the  integral  in  K^x 


To  evaluate  the  integral  in  Eq.  (E44),  we  deform  some  contours  and  hence  attend  to 
branch  points  and  branch  cuts.  As  a  hrst  case  consider  a  simplihed  integral,  namely 

J+  =  (E46) 
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under  the  condition  that  a;  >  0.  The  branch  points  are  ^  =  i  and  ^  =  2i,  both  in  the  upper 
half  .^-plane.  Avoid  crossing  branch  lines  by  running  the  branch  line  up  along  the  positive 
imaginary  axis  of  the  ^-plane.  Let  C  ~  and  ^  —  2i  =  r2e*'^^,  so  we  have 


A _ L  =  !Tp*(A-02). 


Stt 

Y 


<  01,  02  < 


TT 

2’ 


(E47) 


which  implies 

In^— ^  =  ln(ri/r2) +i(0i  -  02).  (E48) 

Deform  the  integration  path  to  run  from  ioo  —  e  down  around  C  =  *  back  up  toward 
ioo  +  e.  Observe  that  in  crossing  the  imaginary  axis  above  Im,^  =  2i,  there  is  no  discontinuity 
in  ln[(^  ~  *)/(C  ~  2i)] ,  so  that  above  ^  =  2i,  the  contribution  from  the  downward  part  of  the 
path  cancels  the  contribution  from  the  upward  part  of  the  path.  Therefore  the  path  can  be 
shrunk  to  a  tight  loop  around  the  line  segment  connecting  the  two  branch  points.  We  have 


(V  1  <  u!  <  2)  In 


C-2t 


In 


^=iw-\-E 


i-2i 


^=iw—e 


0  —  (—271*)  =  271*, 


(E49) 


whence  it  follows  that 


f  271 

J+  =  271*  /  *  dw  =  —  (e-2"  -  e-")  . 

./i  X 


(E50) 


Next  observe  that  the  integral  with  the  branch  points  reflected  about  the  real  axis,  so 
they  are  in  the  lower  half  .^-plane,  is 


JA  =  f  In  =  0,  (E51) 

as  follows  from  pushing  the  integration  path  infinitely  upward.  Thus  Eq.  (E50)  gives  us  also 
the  integral 

/OO  ^2  I  1 

d^  In  =  —  (e-2-  -  e"")  .  (E52) 

e  +  4  a;  ^ 

This  procedure  generalizes  to  produce 

Ji{x)  =  f  ^  _  g-^db)  ,  (E53) 

0-oc  +  1^2  kl 

[Checks  OK  in  MATLAB.] 
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Now  we  evaluate  the  integrals  in  -ftTsi: 


K..(x)  = 


{ki  +  K2  —  Kq){T3ii{x)  —  kQT^i2{x))  +  T3i4(a;)]  + 


kl  +  kl 


T3i3(a;),  (E54) 


where 


T.uix)  = 


In 


e  +  «;i 


=  Ji{x)  =  —  (e 

X 


—  K2X  _  ^  — 


T^12{x)  =  /  di 


00  i^x  p2  I  2 


e^Ki  e  +  «i’ 


T3i3(a;)  =  = 

J_oo 


rp  /  \  aei 

2314(0:)  — 


The  integral  T'3i2(a;)  has  two  terms.  From  the  pole  at  .^  =  mg  comes  a  term 


(E55) 


(E56) 


(E57) 


(E58) 


In  4^  =  i In  4^  +  , 

Kq  Ki  -  Kf^  Kq  \  f^Q-  J 


(E59) 


where  the  second  expression  comes  from  attending  to  the  branch  cut.  The  second  term  of 
7312(2^)  comes  from  a  contour  integral  around  the  line  segment  connecting  the  branch  points 
m2  to  mi  so  that  for  k,2  such  that  branch  crossings  are  no  issue  (that  is,  for  a  <  t\ where 
arg  ^2  =  —ct)  we  have 


i{x)  =  —e  An  ^  ~  27rt  f  d^ 

^0  \  ^0  ~  ^2  /  JiK2 


(E60) 


where 


-2"/ 

JiK2  S  W  f^o  ^0  JiK2 


^  -  mo  c  + 


=  —  {e  (Ei[(fi;i  -  kq)x\  -  Ei[{k2  -  fi;o)a;]) 
_g«:ox  ^  Kq)x]  —  Ei[(fi;2  +  ^0)2^])}  • 


(E61) 


(E62) 


Here  Ei  is  the  exponential  integral: 


e,{z)‘£ 


(E63) 


(For  small  \z\,  Ei{z)  =  —7  —  \n.z  +  z  +  0{z‘^)  where  here  7  is  Euler’s  constant  ~  0.5772 . . . .) 
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As  checked  in  MATLAB,  the  r.h.s  of  Eq.  (E60)  matches  that  of  Eq.  (E56)  for  values  of 
a  less  than  7r/2;  beyond  that,  Eq.  (E60)  gives  the  desired  analytic  continuation  that  would 
come  from  deforming  the  path  in  Eq.  (E56)  to  avoid  crossing  branch  lines. 

The  technique  developed  above  shows  that 


f'lKl 


T3i4(a;)  =  -27r*  / 


J  —IK2 

\  f  K-l  2ki  2\  _ 

—  271  I - \ - ^ — I — ^  )  e 

\X  X-^  J 

27r .  o  o.  ^ . 

— ^ - [Ki  —  K2)  +  U{x)  as  x 

3 


Kn 


—  + 


2K2 


H — 5-  e 


X 


x^ 


:,-K2X 


(E64) 

(E65) 

(E66) 


Putting  this  together  results  in  Eqs.  (6.14)  and  (6.15). 

Note  that  the  correct  branch  of  Ei[{k2  —  Kq)x]  requires  adding  27ii  in  case  the  imaginary 
part  of  K2  —  Hq  is  positive.  This  is  done  in  MATLAB  mlab/branch/K31  .m,  which  exhibits 
smooth  behavior  over  the  range  0  <  a  <  2,  as  shown  by  mlab/pltK3  .m. 


4.  Evaluation  of  Kz2{,x) 


From  Eq.  (E45)  we  have 


^32(2:)  =  - j(G'2  -  Gi), 

Ki  K2 

where  we  define  Rj  by  Ri  ni/  K2  and  R2  *=  ^2/1^1  and 


(E67) 


def 


G.=  /  di 


1  +  «!  -  in^Rj 


+  ^0  +  «!  +  iKiRj 


(E68) 


Then  for  j  =  1,2,  we  split  Gj  into  three  terms: 


Gj  —  Gjo  +  Gji  +  Gj2i 


(E69) 


in  the  following  way.  For  Gj,  j  =  1,2,  deform  the  integration  path  as  shown  in  Fig.  9, 
picking  up  a  residue  contribution  Gjo  at  ^  =  ikq,  where  the  sign  of  kq  is  determined  as 
discussed  in  Appendix  E  7.  Then  split  integration  path  into  two  parts  per  diagram,  one  for 
Gji  around  the  line  segment  in  the  .^-plane  connecting  and  iKj-,  the  other,  Gj2,  is  from 
the  U-shaped  path  vertically  above  Kj. 
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FIG.  9:  Integration  paths  Fji  and  V j2  on  complex  K-plane. 


[To  test  branch  choices  for  the  square  root  and  the  logarithm  for  Gj2-,  observe  that  the 
branch  is  easy  to  hgure  out  for  ^  having  a  large,  positive  imaginary  part.  To  test  the  correct¬ 
ness  of  branch  choices,  we  made  a  scatter  plot  on  the  complex  plane  of  real  and  imaginary 
parts  of  the  square-root  involved  and  of  the  logarithmic  term.  If  both  of  these  scatter  plots 
avoid  branch  lines  and  if  the  branch  choice  is  correct  in  the  limit  of  large  positive  imaginary 
part,  we  are  in  business.  A  neat  thing  about  this  procedure  is  that  it  works  also  when  we 
change  the  integration  variable  for  Gj2  to  y  where  ^  =  i{y  +  Kj).  This  was  done  for  enough 
values  of  arg  K2  and  Kq  to  convince  us  that  the  branches  are  correct  for  the  cases  here.] 
Using  the  technique  discussed  in  Appendix  E5,  we  hnd  Gj2  =  0, 


G,.  =  1 

\  Ka  - 


■  “  tc  ^ 

'4  ^0 


/  kI  —  kI  —  in^Rj 

/k^  —  Kq  +  iK4^Rj 


[and  with 


=  -7rKoe-^°^- 
V 

the  right  choice  of  branch) 


+  «!  -F  in^Rj 


_  1^^ 

fVQ 


—  Hq  —  in^Rj 
kI  —  K,q  +  IK^Rj 


/^Q  —  ^4  —  iK4^Rj 

/kI  —  Kq  -|-  in/^Rj  ’ 

AT7J  o 


(C  +  1)^ 


+  ihaRj 


(E70) 


(E71) 
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In  Eq.  (E70)  the  square  root  has  the  branch  cut  as  the  negative  real  axis,  and  the  square  root 
of  a  positive  number  is  positive.  In  Eq.  (E71)  ^  is  related  to  the  integration  variable  (  by 
^  =  iK,4^{C  +  !)•  The  program  that  calculates  Gj  using  these  equations  is  mlab/Gj  .m.  It  has 
been  tested  against  the  defining  equation  over  the  range  0  <  a  <  1.8,  and  behaves  smoothly 
(no  branch  jumps!)  over  that  range.  That  program  invokes  as  a  subroutine  \kappaOh.m  to 
determine  the  appropriate  value  of  Kq. 


5.  Study  of  a  simplified  integral  relevant  to  Gj  and  to  K2{0) 


Define 


dC 


In 


(E72) 


>-oo  +  1  +  1  +  i 

By  inspection,  X  is  purely  imaginary.  Deform  the  contour  to  come  down  the  left  side  of 

the  imaginary  axis  from  ioo  to  i  and  back  up  the  right  side  of  the  imaginary  axis.  Let 
^  =  i{y  +  1)  to  obtain 


J  = 


dy 


In 


\/l/(l/  +  2)  -  1 


(E73) 


dr  ^/yiy  +  2)  ^/y(y~+^  +  1  ’ 

where  T  runs  from  oo  +  to  the  origin  of  the  |/-plane  and  then  down  around  the  origin 

and  back  out  to  oo  —  iS,  so  the  branch  line  in  the  square  root  is  the  positive  real  axis.  We 
require  that  the  square  root  on  the  lower  side  of  the  branch  line  be  >  0.  From  now  on 
let  \/y{y  +  2)  denote  the  positive  square  root  on  the  lower  side  of  the  branch  line.  Then, 
because  the  change  in  sign  of  the  square  root  from  above  to  below  the  branch  line  cancels 
the  change  in  direction  of  the  path  above  relative  to  the  path  below,  we  have 


J  = 


dy 


In 


Vyiy  +  ‘d) 


'o  \/y{y  +  2)  I  ^yy{y +  2) +  1 


dy 


In 


Vyiy  +  ^) 


In 


+  In 


-Vyiy  +  ^)  ■ 
-\/yiy  +  ^)  + 1, 
Vy(y  +  ‘d)  + 1 


(E74) 


do  ^/y{y  +  2)  y  ^/y{y~+^  +  1  ^y(y  +  2)  -  1  ^ 

Observe  that  \/y{y  +  2)  —  1  has  a  zero  at  y  =  a/2  —  1  which  causes  a  singularity  in  the 
logarithms.  Now  observe  also  that 


In  (^^/y{y  +  2)  -  1  j  =  In  (^1  -  x/y{y  +  2)^  ±  in. 


(E75) 
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Tracing  the  value  of  the  square  root  on  the  chosen  branch  determines  the  sign,  with  the  result 
that  we  find  that  the  path  splits  into  two  pieces  and  that  the  path  beyond  the  singular  point 
dX  y  =  a/2  —  1  makes  no  contribution.  The  plus  and  minus  signs  for  in  turn  out  in  such  a 
way  that  we  obtain 

(in 

I  =  —2ni  /  ,  =  =  cosh~^(v^).  (E76) 

Jo  ^/y{y  +  2) 

This  result  is  confirmed  by  a  MATLAB  integration. 


6.  Integral  for  K2{0) 


The  same  line  of  reasoning  shows  the  relation 


e 


i^a 


In  +  2* 


which  for  the  special  case  of  a  — ^  0  yields 


=  2TTi  I  dye  +  2),  (E77) 

Jo 


+  =  -2iir  f'^dtVW^l 


=  tn 


In  ( 6  +  vT+l^  )  —  bVl  +  b‘^ 


(E78) 


Eq.  (E78)  checks  correct  in  MATLAB. 


7.  Pole  for  residue  contribution  to  Gj 

The  equation  for  the  parameter  kq,  the  location  of  the  pole  in  the  Fourier  transform  of 
the  kernel,  determines  only  Kq.  In  the  defining  equations  for  the  kernel  Kq  enters  only  as 
Kq,  so  the  sign  of  ko  does  not  matter;  however,  in  performing  contour  integrals  a  residue 
from  only  one  of  the  two  poles  in  the  .^-plane  at  ^  =  JzKq  contributes,  so  that  the  sign  of  kq 
matters.  For  the  path  deformations  considered,  the  pole  that  contributes  is  that  specified 
in  the  remark  under  Eq.  (6.9) 
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8.  Evaluating  zeros  in  Fourier  transform  of  electric-field  kernel 


To  evaluate  Eq.  (3.3)  in  cases  of  complex-valued  propagation  constants  and  with  y  set 
equal  to  a,  we  use  the  “kappa’s”  to  put  Eq.  (3.3)  in  the  form 


E(e,a)~-e"Fi(a,e)  +  E2(a,0, 


(E79) 


where,  now  using  “kappa’s”  we  have 


Fi(a,0  i 


«?72  +  ’ 


/OO  -I 

-OO  7i  +  72 

Thus  we  seek  the  values  of  denoted  for  which,  copying  Eq.  (3.3), 


(E80) 

(E81) 


0  —  a)  +  F2(.^,  a). 


(E82) 


where  in  terms  of  “kappa’s”  we  have 

=  2  I  2 1  2  _ 


K. 


2  L 


2  +  2^  ,_2 


Kf  In 


kI  In 


—  2G{Ki,  K,2,  0  +  2G{k,2,  Ki,  ^)] 

+  2Ko  +  2^4^  I  +  0(n^  In u), 


(E83) 


where  u  =  \a\{K  +  |.^|),  with  k  =  max(|fi;i|,  \k2\)  and  we  recall  that  ^4  =  K2\/el{e  +  1)  so 
^2^2 

that  k\  =  ^  ^  ^  ^ ,  and  (7(^2,  ^1,0  is  the  function  defined  in  Eq.  (E41). 

H”  ^2 

Putting  all  this  together  and  dropping  terms  0{u^\nu)  gives 


Fi(^,a)  =  „  ^ 

Ki  -1-  K2 


+ 


2 1  In  -  -  +  ^2  +  ^1)  -  ^1  +  ^2) 

'a  )  nf  —  kI 


IKiK,2K,a 


{kI  -  +  K. 


In  — r^=^= - - In  — r^=^= - — 

/e2^  2  ,  [7^,  2  , 


(E84) 


For  numerical  evaluation,  we  note  that  the  zeros  in  the  .^-plane  are  near  ,^0  = 

±m2\/(e  -f  l)/2  e*'^  where  <5  is  small  and  positive. 
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We  also  find  analytically  (and  checked  numerically): 


•  poo 

F2{ta)  =  -^ - 2/  (E85) 

K-l  K,2  J-oo 

=  -^2^2  (“\/  s'  +  ^'  )  • 

1  2 

The  small-argument  approximation  for  the  modihed  Bessel  function  Ki  (from  [4],  p.  9  and 
[5],  p.  15  for  '^-functions)  is 


2 

.2fc2 


Ki{z)  =  ^  +  ^(ln|+7-^)  +0(^3ln^), 


(E87) 


leading  to,  within  terms  of  order  +  k^)  In  +  k^))  |, 

^2(^,0)  =  -i 


1,1  1_  1  H4  +  -  (^1  + 

o  ^  ^  o 


rvi  rvo 


(E88) 


Eq.  (3.3)  is  evaluated  using  Eqs.  (E88)  and  (E84)  in  the  MATLAB  program  findkOA.m. 


APPENDIX  F:  MATLAB  PROGRAMS 

Taking  inputs  of  e  and  a/h,  the  program  res5.m  calculates  the  hrst  hve  resonant  values 
of  K2h.  For  this  it  calls  resB.m  which  in  turn  calls  nxtB.m.  The  program  resB.m  supplies 
the  error  as  part  of  its  output,  calculated  using  ISD.m.  The  computing  of  the  integral  over 
the  kernel  is  done  by  the  program  ISD.m  and  the  similar  subroutine  ISDC  of  nxtB.m.  These 
call  K.m  which  in  turn  calls  turn  calls  Kl.m,  K2.m,  K3.m.  The  program  K2.m  calls  K20.m 
to  get  1^2(0),  and  K3.m  calls  K31.m  and  Gj  .m.  Many  of  these  programs  call  findkOA.m  to 
obtain  —i  times  the  zero  in  the  .^-plane  of  the  Fourier  transform  of  the  kernel  K{x)]  that 
Fourier  transform  is  accurately  approximated  by  the  subroutine  Ktilde  inside  findkOA.m. 

function [kO] =f indkOA (ep , ell2 , al , a_by_h) 

"/ok  here  stands  for  \kappa 

"/oal  =  -  arg(k)  al=pi/2  for  propagation  constant  real 
"/oCcase  of  most  interest) 

%  refines  function[y]=findk0k2(a,krel) 
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°/o  a  is  |k_2|a  ;  krel  =  k_l/k_2. 
y„gets  zero  from  linear  approx  to  a  complex  function. 

“/owj  =  w(zj) 

zstart  =  i*exp(-i*al) *sqrt ( (1+ep) /2) *ell2* (1+ . 04i) ;  Zrough  starting 
y„approx  to  X=\xi2/ I  k2 1  "2  for  square  of  zero 
del  =  .1; 
zst  =  zstart; 

for  11=  1:7  “/oset  to  1:7 

zl  =  zst+del; 

wl  =  Ktilde(ep,ell2,al,zl,a_by_h) ; 
z2  =  zst-del; 

w2  =  Ktilde(ep,ell2,al,z2,a_by_h) ; 
ykl  =  (w2*zl-wl*z2)/(w2-wl) ; 

zst  =  ykl; 
del  =  del/8; 
end  "/of  or 

testl  =  Ktilde(ep,ell2,al,ykl,a_by_h)  ;yoshows  error  of  found  k_0 
kO=-i*ykl ; 
kOtst=exp(i*al) *k0 ; 
if  imag(kOtst) <0 
k0=-k0; 

end 

function [y] =Ktilde (ep , ell2 , al , xi , a_by_h) 

y„29  July  2009:  fix  KtildeS  of  7  June  2008  for  arb  |k_2| 

y  k  here  stands  for  \kappa 

yal  =  pi/2  for  real  propagation  constant  (case  of  interest) 
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“/oWORKS  FOR  al  near  pi/2;  not  likely  good  otherwise!  I  I 
yoXi=\xi 

y„al  =  -arg(\kappa_2) 

“/oKtilde  used  for  xi  \approx  \xi_0  \approx  .  .  . 
y  i*exp(-i*al)*ell2*sqrt( .5*(ep+l))*exp(i*delta) ,  ... 
y„  delta  small  and  positive. 

X=xi2; 

ph=exp(i*al) ; 
k2=ell2/ph; 
kl=sqrt (ep) *k2 ; 
klsq=kl2 ; 
k2sq=k22 ; 
ftl=kl2+X; 
ft2=k22+X; 

log2=log(X+k2sq)  ;  “/.TEMPORARY  FOR  al  near  pi/2 
logl=log(-X-kl2) -i*pi ; 

gam  =  0.5772156649015328;  %  from  Abramowitz 

“/.function  [y]  =F1  (ep ,  ell2 ,  al ,  xi ,  a_by_h) 
f acl=l . / (klsq+k2sq) ; 

T0=  2*(log(2/a_by_h)-gam) ; 
k4=k2*sqrt (ep/ (ep+1) ) ; 

Tl=(k2sq.*logl-klsq.*log2) . / (klsq-k2sq) ; 
rt3=sqrt (X+k42) ; 

f ac2=i*kl*k2*k4/ ( (klsq-k2sq) *rt3) ; 
lograt2=log( (rt3-i*k2*k4/kl) /(rt3+i*k2*k4/kl) ) ; 
logratl=log( (i*kl*k4/k2-rt3) / (rt3+i*kl*k4/k2) ) -i*pi ; 
T2=f ac2 . * (lograt2-logratl) ; 
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Fl=facl.*(T0+Tl+T2) ; 

“/oend 

y„f  unction  [f  val]  =F2  (ep ,  ell2 ,  al ,  xi ,  a_by_h) 
fl  =  -gam+log(2/a_by_h)+.5; 

F2=  fl+.5*(ft2*log2-ftl.*logl) ./(klsq-k2sq) ; 
“/oond 

y=X.*Fl+F2; 

end 

end 


function  [yg]=Gj (j ,ep,ell2,al,a_by_h,x) 
y„15  July  2009 
y„ell_j=abs  (kappa_j ) 

“/ohere  k  is  for  kappa;  al  =  -arg(k2)=-arg(kl) 
yep  is  dielectric  constant,  ASSUMED  REAL 
k0=f indkOA(ep,ell2,al,a_by_h) ; 
k2=ell2*exp(-i*al) ; 
kl=sqrt (ep) *k2 ; 
if  j==l 

R=kl/k2; 
else  R=k2/kl; 
end 

ell4=ell2*sqrt (ep/ (ep+1)) ; 
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k4=k2*sqrt (ep/ (ep+1) ) ; 
er=10-8; 

rtf=sqrt  (k02-k42)  ;  “/ofrom  study  of  branch  in  12  July  note 
bf  =k4*R; 

G2res  =  i*pi*kO*exp(-kO*x)*log((rtf-bf)/(rtf+bf))/rtf ; 
ulim  =  sqrt (1+R2)-1 ; 

Int_term  =  -2*pi*i*quad(@itg2,0,ulim,er) ; 
yg  =  G2res+Int_terin; 

function  y2=  itg2(t) 

efac  =  exp(-exp(-i*al) * (t+1) *ell4*x) ; 

rrtf =1 . / sqrt (t . * (t+2) ) ; 

tls=(t+l) . "2 ; 

y2=tls . *ef ac . *rrtf . / (tls-exp(2*i*al) *k02/ell42) ; 
end 

end 


function  yf =ISD (ep , ell2 , al , a_by_h, n) 

% compute 

y„\int_02  dx  K(x)  *2*  ( ( (n*pi/2)2-k02)  *  (2-x)  cos  ( .  5*n*pi*x) 
“L-  (4/  (n*pi) )  *  ( (n*pi/2)  2-k02)  *sin(n*pi*x/2) ) 
y„27  July  2009 

k0=f indkOA(ep,ell2,al,a_by_h) ; 
function  yi=itg(t) 

f 1=2* ( ( . 5*n*pi)2+k02) *(2-t) . *cos ( . 5*n*pi*t) ; 
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f2=(4/ (n*pi))*(( .5*pi*n)2-k02)*sin( .5*n*pi*t) ; 


yi=K(ep,ell2,al,a_by_h,t) . *(f l-f2) ; 
end 

er  =  10-9; 

yf=quadv(@itg,0,2,er) ; 
end 


function  [y] =K (ep , ell2 , alpha , a_by_h , xvec) 
T1=K1 (ep , ell2 , alpha , a_by_h , xvec) ; 

T2=K2 (ep , ell2 , alpha , a_by_h , xvec) ; 

T3=K3 (ep , ell2 , alpha , a_by_h , xvec) ; 

y=Tl+T2+T3; 

end 


function [y] =K1 (ep , ell2 , al , a_by_h , x) 

"/X  for  kappa;  21  July  2009 
S=size (x) ; 

rh=sqrt (a_by_h2+x . ~2) ; 
yousmall  =  0; 
k2=ell2*exp(-i*al) ; 
kl=sqrt(ep)*k2; 

yassume  x  a  vector  with  monotonically  increasing  elements 
nv  =S(2);  ynumber  of  elements  in  x 
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if  rh(l)<  .001 


nsmall=min(ceil (nv* . 001/ (10-8+rh(nv) ) ) ,nv) ; 
nbig  =nv-nsmall; 
if  nsmall  ==nv, 
y=Klsmall (rh) ; 
end 

rhsmall  =  rh(l : nsmall) ; 

rhbig=rh(nsmall+l :nv) ;  Zgets  here  only  if  nsmall<nv 
ysmall=Klsmall (rhsmall) ; 
ybig=Klbig(rhbig) ; 
y=  [ysmall  ybig]  ; 
else 

y=Klbig(rh) ; 

end 


function  yf=Klsmall (rrh) 

f ac  =2 . /rrh- (4/3) * (kl2+kl*k2+k22) / (kl+k2) + (kl2+k22) *rrh/2 ; 

yf=pi*fac/(kl2+k22) ; 

end 

function  yf f=Klbig(rrh) 

tl=(l+k2*rrh) . *exp(-k2*rrh) ; 
t2=(l+kl*rrh) . *exp(-kl*rrh) ; 
yff=  4*pi*(tl-t2)  ./((kl4-k24)*rrh.''3)  ; 

end 

end 
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function  [val] =K2 (ep , ell2 , al , a_by_h , x) 
y„accepts  vector  x 
kO=f indkOA(ep,ell2,al,a_by_h) ; 
k2=ell2*exp(-i*al) ; 
kl=sqrt (ep) *k2 ; 

K2null=K20 (ep , ell2 , al , a_by_h) ; 

S  =size (x) ; 
nx=S(2) ; 
er=10-8; 

Tl=x; 

delta=10-7 ; 
for  jj=l:nx 
xtemp=x( j j ) ; 
if  xtemp>delta 

T1 ( j j )=quad(@itg,delta,xtemp, er)+corfun(delta,xtemp) ; 
else 

Tl(j j)=corfun(xtemp,xtemp) ; 

end 

end 

T2=K2null*cosh(kO*x) ; 

y„get  correction  for  starting  integration  from  delta 
function  [cor] =corfun(del ,xc) 
fl=  2*pi/(kl2+k22) ; 
rat=del/ a_by_b ; 
tcorll=log(rat+sqrt (rat2+l) ) ; 
tcorl2=- (2/3) *del* (kl2+k22+kl*k2) / (kl+k2) ; 
tcorl=tcorll+tcorl2; 
tcor2=a_by_h-sqrt (del2+a_by_h2) ; 
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cor=f 1* (tcorl*sinh(k0*xc)+tcor2*cosh(k0*xc) ) ; 
end  “/ocorfun 

val= (kO- . 5* (kl2+k22) /kO) *T1+T2 ; 
function  yi=itg(t) 

yi=sinh(kO* (xtemp-t) ) . *K1 (ep, ell2 ,al , a_by_h,t) ; 

end 

end 


function  yK=K3 (ep , ell2 , alpha , a_by_h , xvec) 
yoell_j=abs  (kappa_j ) 

yaccepts  xvec  as  a  vector  of  values  of  x 
“/ohere  k  is  for  kappa;  al  =  -arg(k2)=-arg(kl) 
yep  is  dielectric  constant,  ASSUMED  REAL 
kO=kappaOh (ep , ell2 , alpha , a_by_h) ; 
k2=ell2*exp(-i*alpha) ; 
kl=sqrt (ep) *k2 ; 
ell4=ell2*sqrt (ep/ (ep+1) ) ; 
k4=k2*sqrt (ep/ (ep+1) ) ; 

S  =size(xvec); 
nx=S (2) ; 
yK=xvec ; 
for  jj=l:nx 

x=xvec ( j j ) ; 

K32=i*kl*k2*k4*(Gj (2)-Gj (l))/(kl4-k24) ; 
KK31=K31 (ep , ell2 , alpha, a_by_h,x) ; 
yK(jj)=K32+KK31; 
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end 


function  [yg]=Gj(j) 
y„15  July  2009 
if  j==l 

R=kl/k2; 
else  R=k2/kl; 
end 

“/olim  =100*pi; 

er=10-8; 

rtf=sqrt  (k02-k42)  ;  y„from  study  of  branch  in  12  July  note 
bf  =k4*R; 

G2res  =  i*pi*kO*exp(-kO*x)*log((rtf-bf)/(rtf+bf))/rtf ; 
ulim  =  sqrt (1+R2) -1 ; 

lnt_term  =  -2*pi*i*quad(@itg2,0,ulim,er) ; 
yg  =  G2res+Int_terin; 

function  y2=  itg2(t) 

efac  =  exp(-exp(-i*alpha)*(t+l)*ell4*x) ; 
rrtf =1 . / sqrt (t . * (t+2) ) ; 
tls=(t+l) . "2 ; 

y2=tls . *ef ac . *rrtf . / (tls-exp(2*i*alpha) *k02/ell42) ; 
end 

end 

end 


function [y]  =K20(ep,ell2,al,a) 
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y„PRIMARY 

“/ok  for  kappa;  a  for  a/h 

y„to  evaluate  K_(x,a)  |  _{x=0}  10  June  2009 

“/ofind  kO  using  =  kappa_0  =  -ik_0(i*kap_l , i*kap2) 

y„which  follows  from  rule  that  kappa_j=-i*k_j 

k2=ell2*exp(-i*al) ; 

kl  =sqrt (ep) *k2 ; 

k0=f indkOA(ep,ell2,al,a) ; 

C=pi* ( l-2*k02/ (kl2+k22) ) / (kl2-k22) ; 
yoTl=2*  (exp  (-kl*a)  -exp  (-k2*a) )  /a ; 
yoT2=2*(kf  (k2)-kf  (kl))/(k0*a)  ; 
yoT3=  kl-k2+lnf  (k2)-lnf  (kl)  ; 
yoK20def =C*  (T1+T2+T3) 
gam  =  0.5772156649015328; 

TlS=k2-kl- (kl2-k22) * (log (a/2) +gam- . 5) /kO ; 

T2S= ( (k22-k02) *log(k2+k0) - (kl2-k02) *log (kl+kO) ) /kO ; 

y=C*(TlS+T2S) ; 

end 


function  [yK31] =K31 (ep , ell2 , al , a_by_h, x) 

“/ocheck;  17July09 

“/ok  for  kappa 

k2=exp(-i*al) *ell2 ; 

kl=sqrt (ep) *k2 ; 

“/oelll=sqrt  (ep)  *ell2 ; 

yyyy 

/o  /o  /o  /o 
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kO=f indkOA(ep,ell2,al,a_by_h) ; 
y„er=10-9; 

yoT311=2*pi*(exp(-k2*x)-exp(-kl*x) )  .  /x; 

y„calculate  T312 

y„if  for  CASE  OF  x  SMALL 

yoa=exp(i*al) ; 

t2=pi* (log( (kl2-k02) / (k02-k22) ) +i*pi) *exp (-kO*x) /kO ; 
if  x<.0001 
el=-log(kl-kO) ; 
e2=-log(k2-k0)  ; 
if  imag(k2-k0) >0 
e2=e2+2*pi*i ; 
end 

elp=-log(kl+kO) ; 
e2p=-log(k2+k0) ; 
else 

el=expint ( (kl-kO) *x) ; 
e2=expint ( (k2-k0) *x) ; 
if  imag(k2-k0) >0 
e2=e2+2*pi*i ; 
end 

elp=expint ( (kl+kO) *x) ; 
e2p=expint ( (k2+k0) *x) ; 
end 

eOx=exp(kO*x) ; 

t5f orin=pi* ( (el-e2) . /e0x-(elp-e2p) . *eOx)/kO;  yas  in  branch/ckT312 
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T312=t2+t5f orm; 

y„  end  T312  (from  Eq.(1.18)  of  20July  res5.tex 
f n2=-k0* (kl2+k22-k02) * (kO*T312/pi) ; 
if  X  <  .001 

tt 1= (2* (kl2+k22) -kl*k2) /3-k02 ; 
tt2=k02-.5*(kl2+k22) ; 
f 1=2* (kl-k2) *tt 1+x* (kl2-k22) *tt2 ; 

yK31=pi* (f 1+f n2) / (kl4-k24) +pi*k0*exp (-k0*x) / (kl2+k22) ; 
else 

fl=(k02-k22+2*kl./x+2./x.''2)  .*exp(-kl*x)  ; 
f2=(k02-kl2+2*k2./x+2./x.''2)  .*exp(-k2*x)  ; 
fnl=2*(fl-f2) ./x; 

yK31=pi* (f nl+f n2) / (kl4-k24) +pi*k0*exp (-k0*x) / (kl2+k22) ; 
yoCHECK  K31  for  values  of  alpha  <  pi/2 
yoK31def =2*quad  (@itg ,  0 ,  pi*70/x ,  er)  /  (kl4-k24) 

"/of  unction  [yi]  =itg(t) 

’/of  acl=cos  (t*x)  .  /(t .  ''2+k02) ; 

y.f  ac2=  (kl2+k22+t .  -2)  .  *log(  (kl2+t .  "'2)  .  /  (k22+t .  ''2) )  -kl2+k22 ; 

yoyi=facl .  *f  ac2; 

end 

end 


function  y=kappa0h (ep , ell2 , alpha , a_by_h) 
“/.ASSUME  0  <=  alpha  <  pi  Modified  18  July  2009 
yoell2  =  abs  (kappa_2)  *h 

y=-i*kh0(ep, exp(i* ( . 5*pi-alpha) ) *ell2, a_by_h) ; 
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yt=exp (i*alpha) *y ; 
if  imag(yt)<0 

y=-y; 

end 

function  khVal  =  khO(ep,k2h,aBYh)yol5  July  2009 

“/ofind  khVal  =  k_2h  where  k_2  is  the  zero  in  zeta-plane  of 

“/oFourier  x-form  of  kernel,  as  approximated  by 

yjai’s  formula  (4.28,11)  of  May  2008,  p24 

“/okhValSq  =  khVal  2; 

“/o  reference  case  is  ep=4. 
rt4=sqrt( .5*(ep2+l)) ; 

E=.5*pi*( . 5*(ep+l)-ep/rt4)/ (ep-1)  . . . 

-i* ( . 5- (ep/ ( (ep-1) *rt4) ) *log( (ep+rt4) / ( l+rt4) ) ) ; 
gamVal  =  -psi(l);  “/oEuler’s  constant  0.5772156  ... 
den=  log(aBYh*sqrt((ep-l)*k2h. ''2/2))+gamVal  -log(2) -i*pi*  . 5 
khValSq=  k2h. ''2 .  *  .  5* (ep+1)  .  * (l-i*E. /den)  ; 
khVal=sqrt_dn(khValSq) ; 

function  y=sqrt_dn(x) 
y=exp(i*pi/4) *sqrt (-i*x) ; 
end 

end 

y„tested  against  k2h=exp(i*phi)  for  phi=linspace(- . 4, 1 . 6) 

“/no  jumps  found,  so  branch  cuts  seem  not  to  cause  trouble. 

I  13  June  2009 


end 


function  kappa2_nxt=nxtB (ep , ell2Sl , alS , inc , a_by_li, n) 

“/oS  for  "start" 
ell2S2=ell2Sl+inc; 
kappa2Sl=ell2Sl*exp(-i*alS) ; 
kappa2S2=ell2S2*exp(-i*alS) ; 
f nl=ISDC (ep , ell2Sl , alS , a_by_h ,n) ; 
f n2=ISDC (ep , ell2S2 , alS , a_by_h , n) ; 
kappa2_nxt=f indX(kappa2Sl ,kappa2S2,fnl ,fn2) ; 
function  z=f indX(zl ,z2,fnl ,fn2) 
z=zl-fnl*(z2-zl)/(fn2-fnl) ; 
end 

function  yf =ISDC (ep , ell2 , al , a_by_h , n) 

% compute 

y„\int_02  dx  K(x)  *2*  ( ( (n*pi/2)2-k02)  *  (2-x)  cos  ( .  5*n*pi*x) 

“L-  (4/  (n*pi) )  *  ( (n*pi/2)  2-k02)  *sin(n*pi*x/2) ) 
y„27  July  2009 

k0=f indkOA(ep,ell2,al,a_by_h) ; 

yokO=kappaOh(ep,ell2,al,a_by_h) ;  Crude  by  reliable  for  test 
function  yi=itg(t) 

f 1=2* ( ( . 5*n*pi)2+k02) * (2-t) . *cos ( . 5*n*pi*t) ; 
f2=(4/ (n*pi))*(( .5*pi*n)2-k02)*sin( .5*n*pi*t) ; 
yi=K(ep,ell2,al,a_by_b,t) . *(f l-f2) ; 
end 

er  =  10-8; 

yf=quadv(@itg,0,2,er) ; 
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end 


end 


num  =  300; 

x=linspace (0,20, num) ; 

KV04=x ; 

KV06=x; 
for  jj=l:num 

KV04(j j)=K(4,l, pi/2, 10-4, x(jj)); 

KV06(j j)=K(4,l, pi/2, 10-6, x(jj)); 
end 

Y04=imag(KV04) ; 

Y06=imag(KV06) ; 
plot(x,Y04,x,Y06,  ’  — 0 
axis([0  20  -.5  .5]) 
xlabel  ( 'k_2xO 

ylabeK' Components  of  kernel  K(k_2x)’) 

title ( ’Comparison  of  kernel  for  a/h=10-4  vs.  10-6’) 


function  [k2hres] =res5 (ep , a_by_h) 
k2hres=ones(l,5)  ;  “/k  is  for  k  not  for  kappa 
err_n=k2hres ; 
for  n=l:5 

val=resB (ep , a_by_h , n) ; 
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k2hres (n)=val (1) ; 
err_n(n)=val (2) ; 

end 

errVec  =  err_n 
end 


function [k2reswE] =resB (ep , a_by_h , n) 
yo28  July  2009  n=n_res 

"/oNEED  more  checking  if  dielectric  constant  ep  complex! 
alS  =  pi/2+.04; 

ell2Sl=.5*n*pi/sqrt((ep+l)/2) ; 
inc= . 02 ; 


for  jj=l:4 

kappa2_nxt=nxtB (ep , ell2Sl , alS , inc , a_by_h ,n) ; 

ell2_nxt  =  abs(kappa2_nxt) ; 

al_nxt  =  argM(kappa2_nxt) ; 

ell2Sl=ell2_nxt; 

alS=al_nxt ; 

inc=inc/7 ; 

end 

kappa2res=kappa2_nxt ; 
k2res  =  i*kappa2res; 

yoerr=0  “/o  speed  up  in  exchange  for  loss  of  error  info 
err=ISD (ep , ell2_nxt , al_nxt , a_by_h, n) ; 
k2reswE= [k2res  err]  ; 
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function  ya  =  argM(z)yogets  al  s.t.  z  =  |  z  |  *exp(-i*al) 


ywith  0  \le  al  <  2*pi 
r=abs (z) ; 
y=imag(z) ; 
x=real (z) ; 
if  y  <=0 

if  X  >0 

ya  =  -asin(y/r) ; 

else 

ya=pi+asin(y/r) ; 

end 

else 

if  X  >0 

ya=2*pi-asin(y/r) ; 

else 

ya  =  pi+asin(y/r) ; 

end 

end 

end  7oargM 

end 
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